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ABSTRACT 


Finite element models can be used to discern the location and severity of damage 
in structures. This is frequently pursued by using the differences in measured and 
computed natural frequencies, in conjunction with the sensitivities calculated using the 
FE model. Given that a modal test produced a limited number of natural frequencies for a 
structure, the concept of Artificial Boundary Conditions (ABC) was developed, which 
yields additional natural frequency information for a structure. This is accomplished by 
artificially imposing additional boundary conditions to the measured data. In this thesis, 
the use of ABC to produce an improved set of structural sensitivities is explored. It is 
shown that the selection of ABC sets is best guided by the strain energy distribution in 


the structure. 
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I. INTRODUCTION 


Finite element models are used to predict the dynamic response of the system to 
various excitations, boundary conditions and parameter changes. When tests are 
performed to validate the FE model, inevitably their results, commonly natural 
frequencies and mode shapes, do not coincide with the expected results from the FE 
model. These discrepancies stem from uncertainties regarding the governing equation of 
the system, simplifying assumptions in the derivation, and inaccurate boundary 
conditions. Clearly one would like to have a better model, based on both the theoretical 
and the experimental results. The problem of how to modify the theoretic model from the 


experimental results is known as model updating (Kenigsbuch and Halevi, 1998). 


Finite Element Models are defined by a large number of physical parameters, but 
a modal test of a structure results in a small number of modal parameters which are used 
to modify the model parameters. That is why updating a finite element model in order to 
produce a reliable, confident prediction of the structural response of a system tends to 


become a difficult task to be accomplished. 


A traditional problem associated with updating a FEM has been obtaining 
sufficient measurements in order to obtain accurate enough results. In theory it should be 
possible to measure as many mode shapes as required. This is not feasible, because the 
available transducers and data acquisition hardware limit the frequency range that can be 
measured. Also, as the frequency increases so does the modal density, causing 
considerable difficulty for the modal extraction algorithms. Even if these difficulties 
could be overcome it would still be impossible to accurately measure the infinite number 
of modes of a structure. A finite element model only accurately predicts the lower third or 
so of the natural frequencies and mode shapes. Therefore the structure could not be 


expected to reproduce all the modes predicted by a model (Ewins, 2000). 


FEM can be very large, extending in many cases to several hundred thousand 
Degrees of Freedom (DOF), or more. Experimental modal analysis rarely uses more than 


a couple of hundred transducers; consequently not all the degrees of freedom in the 
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analytical model will be measured. There are often internal nodes which cannot be 
measured. There are two ways of overcoming this difficulty, namely reducing the 


analytical model order or expanding the measured mode shapes. 


Procedures have been developed for measuring larger experimental databases. 
One method is to identify additional and distinct mode frequencies from the same modal 
test that was already performed to identify the mode frequencies of the standard system, 
without the need for physical modification of the structure. The additional frequencies 
and mode shapes correspond exactly to the mode frequencies found when combinations 
of measured coordinates are constrained to ground. These additional frequencies are 
therefore associated with different boundary conditions for the structure and no physical 


change in the boundary conditions has been made. (Gordis, 1999). 


Since the boundary conditions are not actually physically applied, they are termed 
Artificial Boundary Conditions (ABC). The ABC are imposed on the FEM and 
correspond to ideal constraints. With only one experimental database, this design yields a 
separate FEM for each ABC configuration. With the exception of the imposed boundary 


conditions, each model is identical. 


In a set of spatially incomplete Frequency Response Function (FRF) data, the 
ABC are the constraints that define an Omitted Coordinate System (O-SET). In 
performing a vibration test, a choice is made as to the set of coordinates to instrument 
with response transducers. This set of DOF correspond to the coordinates that are going 


to be actually measured. 


Sensitivity based updating makes use of the modal parameters to help locate a 
difference between a FEM and an existing structure. The parameters experimentally 
determined from the actual structure are the natural frequencies of the structure, the 
corresponding mode shapes, and the modal damping ratios; the parameters needed from 
the initial FE model are the natural frequencies and the associated frequency sensitivities. 
Finite element model sensitivities are the first derivative of the eigen-problem solution set 


with respect to the design variables under consideration. 


Frequency sensitivities are the basis for a linear approximation to compute the 
change in the natural frequencies of a model based on a change in a given design 
variable. Trying to locate the source of this difference, i.e. either the error in the model or 
damage in the structure, is a problem that is generally underdetermined due to the number 
of parameters that can be altered and the limited number of measured parameters. With 
the application of ABC’s a better solution to this problem can often be formulated. This 
method of frequency updating greatly reduces the computational time required to 
calculate the results of a change to large finite element models by eliminating the need to 
resolve the eigen-problem for each iteration. This process has been used extensively for 


model updating (see for example, Flanigan, 1998, Luber, 1995, Dascotte, 1995). 


Furthermore, there are some other approaches to search within the solution 
domain for potential solutions, such as the Constrained/ Unconstrained Optimization 
Routines. The Optimization Routine employed searches for the minimum value of an 
objective function which is constructed to minimize the differences between the 
analytical model parameters and the experimental data. This organized approach of 
searching for the optimal combination of design variables, minimizes the number of 
combinations that are investigated thereby reducing the computational requirements of 


this process. 


The real difficulty to this method of correcting a finite element model is to 
determine which candidate solution most closely matches the experimental data. Multiple 
combinations of design variable changes can result in the same changes to the natural 
frequencies of the model so the accuracy of the frequencies is not the only evaluation 


factor. 


An improved method of solution evaluation is to determine the effects on the 
analytical mode shapes and compare them to the experimental mode shapes. This 
comparison of mode shapes is known as the Modal Assurance Criterion (MAC) and is 
discussed in (Allemang, 1982). The MAC is a measure of how closely the analytical and 
the experimental mode shapes coincide. The MAC has a value from 0 to 1, where the 
candidate solution with a MAC value closer to one when compared to the experimental 


data most closely resembles the experimental case and is selected as the optimal solution. 
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The more meticulous the Objective Function is assembled, that is, the more 
properties of the structure are considered (design variables), the better the error 


predictions are to be found in the Constrained Optimization Routine. 


Three different Objective Functions were evaluated in the present thesis in order 


to correlate the confidence in the damage detection results. 


Il. THEORY 


The matrix equation of motion for a spring mass oscillator (SMO) will represent 


the starting point of the derivation: 
[M] {3} +[C] 3} +[K] 4} = {FO} (2.1) 
Excluding the damping[C], equation (2.1) becomes: 
[M {3} +[K] {x} ={FO} (2.2) 
Rearranging equation (2.2) in a more explicit form, it is written as: 


x, (t) K,, fam K,, x, (t) ff 
a a =i: (2.3) 
M a M,, x, (t) K,, —~ Kn, x, (t) ie 


In 


M, -: M 


where [M] and [K] represent the mass and stiffness matrices which define the physical 


properties of the structure. Both matrices are square and symmetric, that is of sizen xn 


where n is the number of DOF in the model. The vector {F (t)} represents the excitation 


vector of size n x J along the x coordinate. 


Any continuous system has a continuous arbitrary distribution of mass and 
stiffness; hence it also comprises an infinite number of DOF. During modal tests, only a 
finite number of DOF are measured, and these measured coordinates represent the 
analytical set or (ASET), where the measuring instruments (accelerometers) are 
positioned. The remaining unmeasured DOF are considered the omitted coordinate set or 


(OSET). 


A. SPATIALLY COMPLETE AND INCOMPLETE DATA 


It is not physically possible to develop a complete finite element model of a 
structural system over a frequency domain. This would require a number of measuring 
devices (accelerometers) to be equal to the number of DOF of the finite element model, 
hence generating it to have a determined problem where the numbers of equations is 
equal to the number of unknowns, therefore dealing with a set of spatially complete data 
to solve the problem with reliable accuracy. However, practical realities only allow a 
limited number of DOF to be measured in a modal test. Hence the result for a real world 
structure is a measured Frequency Response Function (FRF) matrix that is spatially 
incomplete, involving fewer measured DOF than model DOF, therefore obtaining an 


underdetermined problem. 


For this reason, to perform the structural system analysis, that is either the model 
updating or the damage detection, the analytical mode shapes of the Frequency Response 
Function matrix must be reduced in size such that only the mode shapes identified from 
the test/ measured Frequency Response Function matrix are the ones employed to 


perform the system analysis; this is done via a Reduced Order Model. 


The Reduced Order Model makes use of the ASET and OSET coordinates, where 
as stated before, the OSET coordinates are not associated with the measurement locations 
on the structure, and the reduced Frequency Response Function matrix is composed of 
response information from the limited number of measurements corresponding to the 
ASET coordinates. It is important to note that the reduced FRF matrix implicitly defines 
a dynamically reduced Impedance matrix, furthermore, because the basic system has not 
been altered and that it still has the same number of DOF, even though it was determined 
not to describe all of the DOF, the elements remaining in the reduced FRF matrix are 
identical to the corresponding elements of the full FRF matrix. (Ewins, 1984). The matrix 


still contains all the modal information of a fully described matrix. 


B. ANALYTICAL AND OMITTED COORDINATE SETS 


Rearranging equation (2.3) into the frequency domain for the steady-state 


harmonic response for a full order model without the damping, it becomes: 


ie Oe My, x(t) fi 
ee 0 el ee = (2.4) 
K, ee K., Mi ni Mn xX), (t) ie 


where © is the frequency of the harmonic excitation. Considering that the number of 
DOF’s are comprised for the measured and unmeasured set of coordinates, that is the 


ASET and OSET coordinates, we can express equation (2.4) in the following form: 


Ki K,, 2 M a M Xa tf 
-9 = (2.5) 
K oa Kis M ,, M,, Xo ce 
where the subscript “a” refers to the ASET coordinate and “o” to the OSET coordinate. 


Equation (2.5) is also known as the impedance model and it can be redefined as: 
Leis Loo Xa ea q: 
Lig, Lass) oy. Aa, (2.6) 


where [Z] represents the system impedance matrix evaluated at the forcing frequency, Q. 
Assuming a static constraint, that is no excitation on the omitted coordinates 


then {f,} = {0}, and equation (2.5) is rewritten as: 


Ka K,, 2 M a M Xa ie 
-Q) = (2.7) 
K,, K,, M a M,, xX, 0 
Rearranging equation (2.7) into two separate equations and substituting the OSET 


coordinates in terms of the ASET coordinates, only one equation is required: 


[K.,. Rx} +[K.. hx, }- 27 [LM ., Ke, b+ L., He, j= 0 (2.8) 


Solving equation (2.8) in terms of {x,} it becomes: 
{x,}=[-O?KQM,,]'EKOK,, +O°K 2M. Kx.) (2.9) 
or 


{x,} =[-Z..] 7 [Zou] Oa (2.10) 


By definition the inverse of the OSET impedance matrix[Z,, ] can be expressed as: 


oo 


1 


paar ar. |) = ali cares |e Kolo] (2.11) 


where the Det[e] indicates the determinant and Adj[e] indicates the adjoint matrix. 
From equation (2.11), it is noticeable that the bracketed inverse term does not exist at 


those frequencies which satisfy: 


Det|I -O°K,'M,, |=0 (2.12) 


By definition the eigenvalues which satisfy equation (2.12) are those defined by 


the [K,,] and [M,,] matrices, hence the resulting frequencies correspond to the OSET 
frequencies only. 


Considering that the eigenvalues and eigenvectors of a system are defined by the 
unrestrained DOF of the corresponding system, the OSET coordinates of the OSET 


system are unrestrained and the ASET coordinates are constrained to ground. 


The FRF matrix derivation, as well as the Reduced Order Model are described in 


the following sections. 


C. REDUCED ORDER MODEL 


The process of instrumenting a structure with a finite number of response 
transducers defines a Reduced Order Model, where the impedance of the Reduced Order 
Model is nonlinearly dependent on the impedance of the Full Order Model (Gordis, 
1996). 


In order to perform model updating using spatially incomplete data, a reduction of 
the model data is done, that is the reduction of the analytical model, with no loss of 


generality of the results. 


Considering a full order “exact” FRF model of a structure, we would have a FRF 


matrix of infinite dimensions: 


H H 


oa oo 


H, H 
(l-| = | (2.13) 


where the experimental FRF matrix actually measured in a real test is seen to be a matrix 


partition which has been extracted from the infinite dimension matrix, 1.e. 





la" |= [4] (2.14) 


The overbar in equation (2.14) specifies a reduced model, the superscript x 


denotes an experimentally measured FRF, and the [H,,] matrix represents a structural 


dynamic model that has been reduced using exact dynamic reduction (Gordis, 1996). 


From the partitioned inverse relation of the FRF matrix to its associate impedance, 


ie.[Z]|[H]=1, we have: 
Z Z H H 1 0 
aa oa | aa °| = (2.15) 
Ls Log A, A,, 0 1 
Now, solving for [H/,,| gives: 


ao oo oa 


[Hea] = Ze, —Za.22Z 0) (2.16) 


Since it was assumed no excitation on the omitted coordinates ({f,}= {0} ), the 


exact dynamic relationship between the OSET coordinates and the ASET coordinates 


x, | I 
CEI, | an 


Rearranging equation (2.5) using equation (2.17), the dynamic reduced model 


gives: 


results: 
{f.} =| Zoe ZoPor Zea {Xa} (2.18) 


Comparing equations (2.16) and (2.18), we can see that the partition of the infinite 
dimensional FRF matrix measured in a test is identically equivalent to the inverse of a 
dynamically reduced impedance matrix. Again we see the presence of the OSET system 


dynamics in the term Z,’. The formula of the matrix inverse, given by equation (2.11) 
applies in the same fashion, and since every element in the Z,’ is singular at the natural 


frequencies of the OSET, it is seen from equation (2.16) that the elements of H,' will be 


singular at the OSET natural frequencies (Gordis 1999). 


Given that the measured FRF matrix implicitly defines a dynamically reduced 
impedance model, we pursue the reduction of the FEM in the same manner. That is, we 
calculate an n x n FRF matrix from the FEM, and simply extract the portion of this matrix 
which corresponds to the coordinates measured in a test. Hence the resulting extraction- 


reduced finite element FRF matrix retains all the modal content of the original model. 
D. FREQUENCY RESPONSE FUNCTION [H] 


The Frequency Response Function is a characteristic of the system where the 
frequency response is the frequency domain input-output relationship which can be 
considered the equivalent of the time domain impulse response function (Ewins, 2000). 
The FRF function contains the magnitude and phase (complex amplitude) of the response 


66299 
1 


at the corresponding DOF due to a harmonic force of unit amplitude at DOF “7”. The 
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time domain measured data is transformed to the frequency domain using the Fast 


Fourier Transform (FFT) algorithm. 


The inverse of the FRF matrix [H I. is known as the impedance matrix |Z | , that 


[ZQ)] = [AO] (2.19) 
and as it was implied from equations (2.5) and (2.6): 

[Z(Q)] = | K-2°M + jac] (2.20) 
Analyzing equation (2.6) in terms of an n DOF system 

Zi Zin |[aO) fh 
= tg. a ae (2.21) 
Zn Ue, ies x, (t) ys 
Transferring from physical to modal coordinates let {x}!=[®]{q}, where ® is 


the mass normalized mode shapes of the system and grepresents the generalized 


coordinate set, and rewriting equation (2.21) 
[Z]o}a}= tr} (2.22) 
Pre-multiply by [®]' and expanding [Z] from equation (2.20) 
[] [x][@]-2° [oF [M][o]+ sO] [c][e] ]{q}=[e] {7} 2.23) 


Using orthogonality [o]’ [M Jo]=1, and assuming proportional damping 
[C]=a*[K]+4*[M]: 


[ a7 -2? +2jQe, |{q} =[0] {F} (2.24) 


Where @, is the natural frequency of the i” mode, ¢ is the damping ratio and 
[ @ -? +2 jQa, jis known as the modal impedance matrix and is diagonal. This 
matrix is inverted to find the modal frequency response: 


{x} = lol 


1 
o; -Q? +2jQa, 





[o]' {Ft (2.25) 


Transforming back equation (2.25) into physical coordinates by pre-multiplying 
by [®] and using {3} =[@]’ {F}, the modal decomposition of the FRF matrix in 


physical coordinates is as follows: 


(}-[2]p 


oF -Q +2jQa, | | {F} (2.26) 


where | H(Q)] can be defined as: 


1 T 
[H(Q)]= ela 2)jma | (2.27) 


Writing equation (2.27) in summation form: 





- mod es \o* Vo! y 
Ol py a, -2? +2jQo, een) 


Or for any element in terms of | H i (Q) | : 


_s_fotosy 
[H,(Q)]= 2, o? -9? +2 jo, (2.29) 


E. DRIVING POINT AND TRANSFER POINT FUNCTIONS 


The Frequency Response Function represents a harmonic response in the 
i" coordinate, caused by a single harmonic force applied either at the same i” coordinate 
or at a different j” coordinate. Hence the Driving point Function is one where the 


response coordinate and the excitation coordinate are identical; where as the Transfer 
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Function is one where the response and excitation coordinates are located at different 
DOF points. The FRF is the foundation of the modal testing. It shows a direct connection 
between the modal properties of a system and its response characteristics; hence it 


provides an efficient means of predicting responses (Ewins, 1984). 


The following example is a two DOF system with no damping (Ewins, 2000). It 
displays the characteristics of both a driving point and transfer function for the complete 


FRF. 





Figure 1. | Two-DOF Example. 


Assembling the stiffness matrix: 


For the purpose of this example M,=M,= 1.0 and K, = K,;=0.4 and K,=0.8. 


This yields the following results: 
— 0.7071 — 0.7071 
Mode shapes: {o! \. {o : \ = 
— 0.7071 0.7071 


eo 


Natural frequencies: {a(n ad / sec)} = 0954 
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modes {o! Vo' y 
Applying these values in equation |H ,.(Q)|= ee 
pplying q [77,(Q)] Diao 


For the Driving Point function 


ojo! eile {0.707 1{-0.7074" | {-0.707}{-0.7071" 
2 -Y 04-@? 12-2 





[H,,Q)]= (2.30) 


For the Transfer function. 


outest ites I euteay ile } _ {0.7071}-0.7071}" , £0.7071}{0.7071}" 


i OF -Q? 0.4-0? 1.2-Q? 


(2.31) 
It is interesting to determine what controls whether a particular FRF will have 

positive or negative modal constants and thus whether it will exhibit antiresonances or 

not. As it was noticed from equation (2.29), it is the product of two eigenvectors 

elements, one at the response point and the other at the excitation point. Clearly, if we are 

to consider a driving point, then the modal constant for every mode must be positive, it 

being the square of a number. This means that for a point FRF, there must be an 


antiresonance following every resonance, without exception. This is displayed in the top 


portion of Figure 2 which is the Driving Point FRF [H ‘al : 


The situation for the Transfer Function is somewhat different because clearly the 
modal constant will sometimes be positive and sometimes negative, depending upon 
whether the excitation and response move in phase with each other or not. Thus we 
expect Transfer FRF measurements to show a mixture of antiresonances and minima. 
However, this mixture can be anticipated to some extent because it can be shown that, in 
general, the further apart are the two points in question, the more likely are the two 
eigenvector elements to alternate in sign as one progresses through the modes (Ewins, 
2000).This is shown in the lower plot of Figure 2 which is the Hj2. Conversely in the 
region between the two natural frequencies the two terms are additive and thus do not 


create an antiresonance. 
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Figure 2. | Two-DOF Frequency Response Function. 


The principles demonstrated in the above example can be applied to any number 
of DOF. This concept of the existence of an antiresonance between any two modes of a 
driving point FRF is of great significance in the use of Artificial Boundary Conditions as 


it is shown in the next chapter. 
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HI. ARTIFICIAL BOUNDARY CONDITIONS 
CONFIGURATION OF NATURAL FREQUENCIES 


The improvement of a FEM is often a necessary step in order for the model to be 
used with confidence in the prediction of structural response. This improvement is 
difficult to achieve due to the spatially incompleteness of the measured FRF matrix 
giving an underdetermined problem. Undoubtedly there exists a recognized need to 
increase the size of the known parameter data base; it has been found that a number of 
additional mode frequencies can be identified from the same modal test without the need 
for any physical modification of the structure. These additional and distinct mode 
frequencies correspond exactly to those mode frequencies found when combinations of 
measured coordinates are restrained to ground. This last idea is referred as the ABC 
concept, which as it will be shown, improves the updating process by associating the 


natural frequencies with different boundary conditions. 


A. RELATION BETWEEN THE DRIVING POINT ANTIRESONANCES TO 
ABC FRQUENCIES 


Evoking the Driving Point Function concept discussed in the last chapter it can be 
shown that the ABC’s frequencies of a given ASET coordinate are defined by the 
corresponding driving point antiresonances. Furthermore, the ABC’s mode frequencies 
not only provide a greater number of frequencies for the system, but also provide a means 


to reduce ill-conditioning in the solution of the sensitivity equations. 
From equation (2.29) the driving point FRF is 


modes (@')’ 
rl @ —O? 





H,(Q)= >) (3.1) 


where ®; is the mass normalized mode shape element, w, is the kth natural frequency, 


and Q is the forcing frequency. 
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The frequency of the anti-resonance of H,,(Q) is given by 


12 2.2 
2 _ Re, +R, 


anti—res Ri, +R), (3.2) 


where the modal residue is giving by Ri‘ = {M(:, k)} {@(:, k)}. The following example 
(Gordis, 1999), shows how frequency is identical to the natural frequency of the system 


in Figure 3, to the driving point DOF constrained to ground as shown in Figure 4. 


1. Calculation of ABC Frequencies for a 2-DOF System 


Given the shown undamped system of Figure 3: 


k1 k2 


Figure 3. | Two-DOF system. 


Let M; = M2= 1.0 and K; = K2= 1.0, the frequency of the single antiresonance of 


the system isQ = 2 rad/sec, which is identically equal to the single natural 


antires 


frequency of the ABC system of Figure 4 which @= V2 rad/sec. 





Figure 4. | ASET, DOF #1 constrained to ground. 
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Figure 5. — Driving Point FRF, Hy:. 


The H,,(©) plot in Figure 5 graphically shows the calculated mode frequencies, 
as well as the anti-resonance frequency as it was computed from equation (3.2) is 


QO = +2 rad/sec. or 1.41 rad/sec. 


antires 


By the use of equation (2.16), the frequencies of the ABC system can be found 
just by calculating inverse of H,,(Q), where the ASET is pinned at DOF #1, and H,'(Q) 
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is the inverse of H,,(Q). A clear correlation is depicted in the bottom graph of Figure 6, 


where the single natural frequency corresponds to the antiresonance frequency. 


Furthermore, this singular frequency happens to be the natural frequency of the ABC 


system obtained by constraining DOF #1 to ground (Figure 4). Hence it is shown the 


power and usefulness of equation (2.16), where the driving point antiresonance 


correspond directly to the natural frequencies of the system, with the driving point 


coordinate constrained to ground (Gordis, 1999). 


Plot A. Driving Point FRF, H11 
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Figure 6. Plot A. Driving Point H11(Q) of system 1, Plot B. [Ha(Q)]' of system 2. 
2. Calculation of ABC Frequencies for a Free-Free Beam 


To further show the usefulness of the ABC concept, a free-free beam is analyzed 


(Gordis, 1999). The present model depicted in Figure 7 is examined, it consists of 10 


beam elements and each element contains two nodes with two DOF each, giving a total 


of 22 DOF system. The translational DOF correspond to the odd numbers where as the 
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rotational DOF correspond to the even numbers. Response transducers are set up at the 


translational DOF: 1, 5, 9, 13, 17, and 21 as follows: 





Figure 7. Transducers located at DOF 1, 5, 9, 13, 17, and 21. 


The set of coordinates chosen to be measured with response accelerators represent 
the ASET, defined as [1 5 9 13 17 21]. It can be assumed that the excitation is applied at 
each of the ASET DOF, resulting in a 6 x 6 FRF matrix. The impedance matrix H;'(Q), 
the scalar inverse of H,,(Q), is calculated at excitation frequencies from 0-800 hertz. 


Figure 8 shows the driving point FRF of the system. 








7 
ie} 100 200 300 400 500 600 700 800 
Frequency (Hz) 











Figure 8. Driving Point FRF ASET DOF: 1 5 9 13 17 21. [From Gordis, 1999]. 
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The ABC system frequencies are found by calculating the impedance matrix 
H,'}(Q). By plotting the 1,1 element H,,(Q) of the matrix, the ABC system frequencies 


can be seen in Figure 9: 








Frequency (Hz) 








Figure 9. Impedance FRF for ASET DOF: 1 5 9 13 17 21. [From Gordis, 1999]. 


The ABC frequencies correspond exactly to the natural frequencies of the system 
with all of the measured coordinates ASET constrained to ground. This configuration is 


shown in Figure 10: 





Figure 10. ABC Configuration ASET [1 5 9 13 17 21] (Measured coordinates 
restrained to ground). 
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IV. SENSITIVITY BASED FEM UPDATING USING ABC 


The disparity in the number of known parameters (measured) versus the number 
of parameters to be adjusted in order to update a FEM, defines a spatially incomplete 


structure as well as an underdetermined problem. 


The importance of obtaining sensitivities of the dynamic characteristics of a 
system with respect to the changes of its parameters lies in the fact that they are critical 
for achieving efficient design modification, and for predicting the optimal structural 
modifications within a prescribed dynamic characteristics change. In such a case, there 
will be numerous possible modifications that can accomplish the change. Sensitivity 


analysis can provide the most effective answer. 


From a mathematical point of view, sensitivity analysis implies differentiation. It 
can be carried out on the modal data from the FRF of a dynamic system. The variables in 
this analysis are usually system parameters such as mass or stiffness parameters. The 
functions of the analysis can be natural frequencies and mode shapes of the system 
(modal data), or its FRF. The assumption of infinitesimal changes demanded by 
differentiation applies to sensitivity analysis. Hence the sensitivity based updating 
concept is used for error localization in order to find the parameters requiring 
adjustments. The ABC concept can be used in addition to the standard baseline system 
mode frequencies in model updating and damage detection. As stated previously, the 
ABC frequencies correspond to the same structural system as the baseline frequencies, 
but with different boundary conditions. The governing equation for sensitivity based 


updating is: 


{Aw} =[T]{ADv} (4.1) 


where {Aa} is a vector of natural frequency errors. The errors are the difference 


between the experimental and the analytical natural frequencies {@* —w*}. The {ADv} 
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term is the vector of changes to be calculated for the specified model parameters, known 


as the design variables, and [7] is the sensitivity matrix. 
Each column of the sensitivity matrix represents an element of the model, where 
as each row represents a mode of the model. 


Each ABC system defines additional rows to equation (4.1) that is it adds more 


modes, defining the equation: 


{Aa [7 {abv} (4.2) 


da; 
where 7° =—+— and @ 
ADV, 


represents the ith natural frequency of the kth ABC 


configuration system. The baseline system quantities are represented with the superscript 
“0” and the ABC systems are identified with the superscripts from | to “k”’, where k 


represents the number of ABC used in the system. Combining baseline and ABC system 


equations: 
2(0) (0) 
AQ@pase T BASE dV. 
Ao TO 1 
; ABC1 = : ABC1 : (4.3) 
f : : : dV, 
Aq?! ) T! ) 


ABCk ABCk 


The degree of coupling between the ABC systems and the baseline system in 
equation (4.3) can be adjusted by deleting or retaining individual columns of the ka |. 


Partial coupling between the baseline system and the ABC system is established if the 


design variables (DV) are partitioned. That is, some of the DV are associated only with 
the baseline system, {ADV°}, some are associated with both the baseline and the ABC 
system, {ADV°":, and some are associated only with the ABC system, {ADV'}, (Gordis, 
1999).The use of the ABC system sensitivities helps to eliminate or at least to reduce the 


problem of poorly conditioned or rank-deficient [T] matrices. Columns of Fag can be 


replaced with columns of [r* | in order to improve the conditioning. Two closely spaced 
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elements in a model will have nearly dependent columns of Ear preventing the 
discrimination of error or damage between two elements. One column of [re associated 


with one of the two elements can be replaced with a column of [r' ] , and the associated 
baseline system natural frequency replaced with the ABC system frequency (Gordis, 
1999). 


A. SENSITIVITY MATRIX DERIVATION [T] 


Several approaches are available for performing a sensitivity analysis. The 
following derivation is based on the parameterization of the eigenvalue problem (Inman, 


1989). Consider an n DOF system defined by: 
M(DV)x(t)+ K(DV)x(t) =0 (4.4) 
where the corresponding eigenvalues problem is: 
[K —AiM ]{@;} = {0} (4.5) 
Differentiating equation (4.5) with respect to the design variable {DV}, as it is the 
vector of design parameters representing the change in mass matrix [M | and/or change 


in stiffness matrix| K | , which are adjustable for the FE model. 








AK AM Adi AQVi 
i ar Se top HK -amy] 2° = (9 (4.6) 


Expanding equation (4.6) and pre-multiplying by {py leads to: 





(0) | fo —afoy] jo} | fey tate) 0)" tx-any S| 10) 4.9 


ADV ADV ADV 


And invoking the matrix identity property {a\' [b] {ch = fel" [b]{a} the last term 


on the left hand side can be replaced by: 





a | [K -AiM]{@:\ =0 (4.8) 
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Since [K —4,M ]{@,} = {0} the overall equation is reduced 





(0) | AE [toy -afoy"| SE |to}-| 2 |foy" to} (0) 9) 


Using orthogonality, where [@]’ [M |[®:]=1 equation (4.9) further reduces to 





(0)"| A |toy-afooy"| SE ]oj-[ 22 [- 10 


Bringing the mass and stiffness portion to the right hand side 


Adi |-cor"| AK _, AM |e (4.11) 


ADV ADV ADV 








Finally it is demonstrated that the sensitivity matrix involves differentiation with 
respect to the design parameter where as for this case are either consider mass and/or 
stiffness, and the function of the analysis is with respect to the natural frequencies of the 


system as it is shown in equations (4.12) and (4.14): 








Ae |-(09"| SE [to where [AK] =[Kx- Kz] (4.12) 

and 
[Tigres | =(0)"] lt (4.13) 
er | 03"|-a at [toy where [AM ]=[Ms— Mb] (4.14) 

and 
au] ={04" | 2-H [f} 4.19 


B. BEAM SPECIFICATIONS 


A Finite Element computer program was written to calculate the Sensitivity 
matrix of a twenty element cantilever beam model, which consists of an aluminum bar, 
with a partially drilled and tapped 10/32” diameter hole at the free end, which allows the 


load cell to be attached to the structure. Table 1 provides a summary of the beam 




















specifications. 
PARAMETER VALUE 
Length 42 inches 
Width 1.5 inches 
Thickness 0.5 inches 
Density 0.11 lb/ in? 
Elasticity Modulus 10 E6lbf /sec’—in. 














Table 1. | Beam Specifications. 


C. MASS AND STIFFNESS SENSITIVITY ANALYSIS 


The computer program built to assemble the sensitivity matrix calculates the 
mode frequency sensitivities according to equation (4.11). This process is performed in 
two steps. The first step calculates the Mass Sensitivity matrix and the second calculates 
the Stiffness Sensitivity matrix. Hence the [T] matrix is of size m x p, m being the 


number of modes and p being the number of elements. 


In the Mass Sensitivity calculation, a small perturbation of 1% of mass was 
applied to each element of the beam; this is done using equation (4.14). Each subsequent 


column of [T] was calculated at every element of the model. 
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On step two for the Stiffness Sensitivity matrix calculation, the same 1% stiffness 
perturbation is applied iteratively to each element and equation (4.12) is used to assemble 


the corresponding part of the Sensitivity matrix. 


Once the two steps are completed, both sensitivities matrices are assembled 
[T(M)] & [T(K)], hence obtaining a [T] total of size m x 2p, where the first set of twenty 
columns yield the changes in mass and the second set o twenty columns yields the 


changes in stiffness. 


The following plots show a normalized mass and stiffness sensitivity distribution 
along the 20 elements of the beam model. Each element is indicated by a bar which 
represents its corresponding response to the change of the system natural frequency, with 


respect to the corresponding change of the design variable. 


It will be seen in an overall sense on all of the Figures below how significantly the 
sensitivity is influenced by the mode shapes. It is clear from Figure 11, how the Base 
Stiffness Sensitivity is more likely to detect a change in the natural frequency closer to 
the fixed end of the cantilever beam model, rather than at the free end, where the red 


sensitivity bars are almost zero. 
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Figure 11. Stiffness Base Sensitivity Matrix. Modes | to 5. 


In Figure 12, it is evident how the Base Mass Sensitivity is more likely to detect a 
change in the natural frequency, closer to the free end of the cantilever beam model, 


rather than at the fixed end, where the magenta sensitivity bars are almost zero. 
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Figure 12. Mass Base Sensitivity Matrix. Modes | to 5. 


It should be noted that the Mass Sensitivity rows are represented in absolute 
values, therefore an increase of mass would generally lower the system natural 
frequencies, but, it was preferred to show the influence with respect to element position 


and hence all magnitudes were made positive. 


The next four figures show a particular trend followed by the Stiffness and Mass 


Sensitivities when ABC’s are applied to the system. 
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The green pointer with triangular shape at element twelve of Figures 13 and 15 
and element four of Figures 14 and 16, represent the pinned ABC applied at DOF 25 and 


9, respectively as follows: 
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Figure 13. Stiffness Sensitivity distribution with ABC at DOF 25. 
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Figure 14. Stiffness Sensitivity distribution with ABC at DOF 9. 
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Figure 15. Mass Sensitivity distribution with ABC at DOF 25. 
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Figure 16. Mass Sensitivity distribution with ABC at DOF 9. 
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In general when ABC’s are applied to the Stiffness Sensitivity matrix, the model 
has a tendency to increase the values of their sensitivities at that particular “pinned” 
ABC, whereas for the Mass Sensitivity, the pinned ABC seems to drive the sensitivity 
values away from the DOF’s where fixed conditions are emulated; these fixed conditions 
can be thought as those representing the fixed end of the beam as well as the pinned ABC 
set up at the corresponding DOF’s 25 and 9 as it was showed on Figures 15 and 16, 


respectively. 


In summary, the previous Figures expressed that in places of low sensitivity, the 
beam model is not very likely to predict the changes in natural frequencies due to the 
change in the design variable, and that is, in a real world scenario those low sensitivity 
areas would be the places where damage, if any, is not going to be accurately detected. 
Hence, the sensitivity concept leads to a key role in the updating of a FEM. Furthermore, 
the sensitivity analysis showed us two important points. The first is the big influence of 
the mode shapes into the sensitivity distribution of the model, and two; how the correct 
manipulation of the ABC’s especially with regard to the stiffness sensitivity portion, 


could lead to a better prediction of the errors in the FEM updating. 
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V. KINETIC AND STRAIN ENERGY CORRELATION TO 
SENSITIVITY MATRIX 


Because the natural frequencies and mode shapes of a structure are dependent on 
the mass and stiffness distributions, any subsequent changes in them should, 
theoretically, be reflected in changes in the frequency and mode shapes of the structure. 
Hence, the first step in a vibration problem is to determine the important natural 
frequencies of the system; there are two general methods to determine the natural 


frequencies of the system: The Newtonian and the Energy method approach. 


The Rayleigh’s Energy method employs an energy balance which is essentially 
the scope of this chapter. For the sake of the following development, the Euler-Bernoulli 


beam theory is utilized, where the Euler-Bernoulli equation for beam bending is: 


ov O° Ov 
EI = q(x,t 5.1 
Q ae ae | ae q(x,t) (5.1) 





where v(x,f) is the transverse displacement of the beam, @g is the mass density per 
length, EI is the beam rigidity, g(x,t) is the externally applied pressure loading, t and 
xindicate time and spatial axis along the beam axis. Rearranging equation (5.1) to 
develop the finite element formulation and the corresponding matrix equations, the 


averaged weighted residual of equation (5.1) is: 


i=) poe. eyo” )_ odx = 0 (5.2) 
Paz ael” ae) 4 


where Lis the length of the beam and @ is a test function. The weak formulation of 
equation (5.2) is obtained from integration by parts twice for the second term of the 


equation. In addition, discretization of the beam into a number of finite elements gives: 





px S| foStours fart va a -f aoa s(-Ver use) =0 (5.3) 


2. 
i=] \ Qe Ox 
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where the slope is represented by: 
0= [=| (5.4) 


Bending moment: 


Orv 
M =EI (= (5.5) 
Shear force: 
V =-El (S (5.6) 
and the Loading: 
w= (= (5.7) 


Q* represents the element domain and n is the number of elements for the beam. 


Shape functions for the spatial interpolation of the transverse deflection v in 
terms of nodal variables are considered. Elements are measured as having two nodes, one 
at each end as shown in Figure 17. The deformation of a beam must have continuous 
slope as well as continuous deflection at any two neighboring beam elements. To satisfy 
this continuity requirement each node has both deflection, v, and slope @ as nodal 
variables. In this case any two neighboring beam elements have common deflection and 
slope at the shared nodal point. This satisfies the continuity of both deflection and slope. 
The Euler-Bernoulli beam equation is based on the assumption that the plane normal to 


the neutral axis before deformation remains normal to the neutral axis after deformation. 


Ov). pines fos 
This assumption denotes 0 = [=| (i.e. slope is the first derivative of deflection in terms 
x 


of x ). (Kwon and Bang, 2000). 
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Al a2 





X1=0 X2=L 


Figure 17. | Two-noded beam element. [After Kwon and Bang, 2000]. 


Because there are four nodal variables for the beam element, we assume a cubic 


polynomial function for v(x): 
V(X) =C, +O, X46,X° +6,x° (5.8) 
The slope is computed from equation (5.8): 
O(x) =c, +20,x+3¢,x" (5.9) 
Evaluation of the deflections and slopes at the nodes yields: 


v(0) =c, =, 
A(0) =c, =8, 
vl) =c,+¢l+o0 +oP =v, 


O(1) =c, +2c,1 +3c,l’ = 8, 


(5.10) 


where / represents the element length. Solving equation (5.10) for c, in terms of the 


nodal variables v, and 6@,, substituting the results into equation (5.8) gives: 
V(x, 0) = Hy (x), (0) + A, (x)6, (0) + A, (xv, (1) + H,(9) 6, (0) (5.11) 
where H,(x) represents the Hermitian shape functions, which are of C' type, meaning 


dv . : ; 
that both, v and a are continuous between the neighboring elements. 
x 
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In general, C” type continuity means the shape functions have continuity up to 
the n” order derivative between two neighboring elements. The Hermitian Shape 


functions for the two node beam element depicted in Figure 17 are: 


Ho) =1 3 
H)=1- 745 
3x 2x8 on) 
Rae Ege 
H,(x)= Haz 








0 


HL &% 





Figure 18. Hermitian beam element. [From Kwon and Bang, 2000]. 


A. STRAIN AND KINETIC ENERGY DERIVATION 


Regarding the Strain Energy due to bending, consider a differential element of 


beam of length dx , where the axial stress due to bending is: 


ga (5.13) 


38 


where / is the cross section area of inertia, the Strain Energy per unit volume is: 
1 
aa a (5.14) 


where di is the change in length of a fiber of the beam segment dx. The work done by 


the applied moment induces axial strain in the fiber, which is stored as potential (strain 


energy). 


Let F =odA, and dL= SL = €dx, then the Srain Energy becomes: 


dU = +(odA) (ed) = 5 oedAds = 5 oedVol (5.15) 
Let dVol =1, then equation (5.15) becomes: 


dU = + osdVol =—08 (5.16) 


: M: : : , 
Since e=T and from equation (5.13), the Strain Energy per unit volume 


becomes: 
1 (Mz) 
aw -3(%) (5.17) 
2E\ 1 
The total Strain Energy is: 
L 2 
1 (M 
U = fav =| f( | dAdx (5.18) 
vol OA 2E I 


Using equation (5.5) and J= } z dA, the Strain Energy stored in a beam of 
A 


uniform cross section that has a deflected elastic curve v(x) is: 


_EIF 0’ v(x,f) . 
ea i a dx (5.19) 


0 
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With respect to the Kinetic Energy, let the displacement for any cross section 


along the beam, to vibrate in a single mode at its corresponding natural frequency a, : 
v(x,t) = H,(x)*C, sin(a,t + @,) (5.20) 


Where C,sin(@,t+6,) represents the temporal nodal coordinates comprised of 
displacements v,(t) or slopes 6,(t) as it is implied from equation (5.11). The velocity of 


the cross section becomes: 
v(x, t) = H,(x)* @,C, cos(a@,t+ @,) (5.21) 
Let C, =1, then the maximum velocity is: 


(x)= A,(x)*@ (5.22) 


Vmax 


Notice that v(x) becomes a function of x only. And the Total Kinetic Energy 


max 


for the beam is (mass/length y ): 


2 


L 
T,., =~mv2,, =~ [(H,(x)*a,) dx 5.23 
max Vinax 71 ( \ )*o, ( ) 


Equations (5.19) and (5.23) were used to derive the Strain and Kinetic Energy 
Sensitivities, respectively in the twenty element beam which characteristics are described 
on Table 1. In order to derive the Strain Energies, in accordance with equation (5.19), the 
second derivatives of the Hermitian shape functions were multiplied by the 


corresponding displacement v, or rotation 0. of the eigenvector matrix (mode shapes), 


generated from the Build2beams.m program in MATLAB. Once the iterative process is 
completed for the number of elements, the Strain Energy matrix is assembled and it is of 
size m Xx p, where the m rows represent the displacement and rotations and the p columns 
represent the different mode shapes, as opposed to the Sensitivity matrix where the m 
rows represent the mode shapes and the p columns represent the different elements of the 


model. 


In regard to the Kinetic Energy calculation, it is implied from equation (5.23), no 


modifications were performed to the Hermitian shape functions which are multiplied by 
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the natural frequencies generated in the Build2beams.m MATLAB program. After the 
iterative process is completed, the Kinetic Energy matrix is assembled with the same 
characteristics of the Strain Energy matrix, that is the m rows representing the 


displacements and/or rotations and the p columns representing the different mode shapes. 


The following graphs show the normalized plots for Strain and Kinetic Energies 
distribution, in the same style as the Sensitivity distribution figures shown in the previous 


chapter. 


Again each subplot corresponds to a specific mode shape and each bar represents 
the level of either Strain or Kinetic Energy stored at each element along the 
corresponding mode shape. The red bars on Figure 19 depict the Strain energy 
distribution along the twenty elements Base system of the beam over the first five mode 
shapes, on this graph we can see how the Strain Energy exactly correlates to the Stiffness 


Sensitivity of the cantilever beam model, as it was shown in Figure 11, Chapter IV. 
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Figure 19. Strain Energy distribution base system. 


The magenta color bars in Figure 20 depict the Kinetic Energy distribution along 
the twenty elements Base system of the beam over the first five mode shapes, on this 
graph we can conceivably see how the Kinetic Energy is exactly correlated to the Mass 


Sensitivity of the cantilever beam model, as it was shown in Figure 12, Chapter IV. 
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Figure 20. Kinetic energy distribution base system. 


In the same way, a comparison of the Strain and Kinetic Energies to the Mass and 


Stiffness Sensitivities is performed with the ABC’s applied. 
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Figure 21 depicts the Strain Energy distribution along the twenty elements with 
ABC applied at DOF 25 of the beam over the first five mode shapes, in this graph we can 
see how the Strain Energy distribution is exactly correlated to the Stiffness Sensitivity 
distribution of the cantilever beam model with the ABC applied at the same DOF, as it 
was showed in Figure 13, Chapter IV. 
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Figure 21. Strain energy distribution with ABC at DOF 25. 
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Figure 22 depicts the Kinetic Energy distribution along the twenty elements with 
ABC applied at DOF 25 of the beam over the first five mode shapes, on this graph we 
can see how the Kinetic Energy distribution is exactly correlated to the Mass Sensitivity 
of the cantilever beam model with the ABC applied at the same DOF, as it was showed in 
Figure 15, Chapter IV. 
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Figure 22. Kinetic energy distribution with ABC at DOF 25. 
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It has been proven for the cantilever beam modeled, that the Strain and Kinetic 
Energy matrices are directly correlated to the Stiffness and Mass Sensitivity matrices, 
respectively; therefore we can conclude that the same trend observed on the Stiffness 
Sensitivities distribution from the last chapter applies for the Strain Energy distribution, 
as well as the same trend observed for the Mass sensitivity distribution applies for the 


Kinetic Energy distribution. 


It is also seen how the Stiffness Sensitivity matrix is directly correlated to the 
flexural properties of the structure and how the Mass Sensitivity matrix is directly 
correlated to the Kinetic energy of the structure, which in other words indicates it is a 
function of the velocity stored in the structure at each specific mode shape. That is why 
the Figures depicted of Mass Sensitivity and Kinetic Energy distributions showed low 
sensitivity at the hard points of the structure (i.e., high stiffness areas) and higher 
sensitivity on the portions away from the fixed points, where there is minimum Potential 


energy, but maximum Kinetic energy. 


This is going to play an important role as it will be seen in the next chapter, as it is 


a very useful tool in the endeavor of predicting damage detection in the model. 
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VI. DAMAGE DETECTION 


It was shown that there exists a direct correlation between the Mass sensitivity 
with respect to the Kinetic Energy matrices as well as the Stiffness sensitivity with regard 
to the Strain Energy matrices. Furthermore with the tendencies observed and discussed at 
the end of the last chapter when the ABC’s were applied; several scenarios were run in 
order to find the best choice of ABC coordinates, such that it could result with the best 


error estimate of the structural model. 


As the number of elements is augmented, the finite element model accuracy 
improves. The disadvantage is that as the number of elements increases, so does the 
computational time. The twenty element model was used to perform the error prediction 
calculation in the present thesis, making use of the first five modes to perform the 
analysis, as it was decided due to the reliable convergence, when the number of elements 


of the model gets increased, as it can be seen on Table 2. 









































# OF ELEMENTS 10 20 42 
MODE | 8.632305 8.632298 8.632298 
MODE 2 54.09948 54.0978 54.09769 
MODE 3 151.5137 151.4776 151.4752 
MODE 4 297.1136 296.8493 296.8317 
MODE 5 491.9194 490.7656 490.6868 
MODE 6 736.9511 733.2692 733.0091 
MODE 7 1033.978 1024.51 1023.809 
MODE 8 1385.246 1364.733 1363.099 
MODE 9 1791.106 1754.315 1750.902 
MODE 10 2226.603 2193.797 2187.248 

















Table 2. Frequency (Hz) versus beam elements. 


In order to run the experiments, two FEM beams were modeled in MATLAB, one 
being the analytical model and the other being the experimental, both models are exactly 
the same and correspond to the cantilever beam depicted in Figure 23, the beam 


specifications are those stated on Table 1. 
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The experimental beam was used to introduce a known perturbation at a known 


element as it is requested from the MATLAB computer program. 


The errors implemented into the experimental beam consisted of a 10 percent of 


either mass or stiffness perturbations. 





20 Element Cantilever Beam. 


Figure 23. 


A. EXPERIMENT 


In order to find a set of ABC’s, such that we could produce a set of modes with an 
improved sensitivity distribution, consequently a better error estimate; twenty four 
different cases were run to visualize the error perturbations introduced on the 
experimental cantilever beam model. Twelve scenarios consider the analysis of one ABC 
and one damaged element, and twelve scenarios consider the analysis of multiple ABC’s 


and multiple damaged elements. Table 3 illustrates the experiment set up as follows: 
































ONE ABC/ ONE DAMAGED MULTIPLE ABC’S/ MULTIPLE 
ELEMENT DAMAGED ELEMENTS 
TYPE OF DAMAGE NUMBER OF CASES TYPE OF DAMAGE NUMBER OF CASES 
STIFFNESS 6 STIFFNESS 6 
MASS 6 MASS rs 
Table 3. | Twenty four Experiment Cases. 


On the following cantilever beam drawings, the pinned ABC’s are represented by 
the red triangles, the elements with the stiffness perturbations are represented in brown 
color and the elements with the corresponding mass perturbations are represented in 


magenta. 
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The selection of the ABC’s was chosen following the sensitivity patterns 
identified and discussed on the last two chapters in relation to the mass and stiffness 


behaviors with respect to the Kinetic and Strain Energies performance, respectively. 
The one ABC/ one damaged element cases are analyzed with 3 different graphs: 


e The cantilever beam with the graphic error representation of the error, as 


well as the pinned DOF where the ABC is applied, 


e The stiffness or mass sensitivity distribution along the beam for the 


specific case in analysis, and 


e The actual damage detection (error prediction) of the model. 


B. STIFFNESS DAMAGE DETECTION WITH ONE ABC 


iB Case 1. Stiffness Error at Element 5, and ABC Pinned at DOF 9 






— (Sc ER SON La NR Dem em 





Figure 24. Stiffness error at element 5, ABC pinned at DOF 9. 


Figure 25 shows the Stiffness sensitivity distribution along the beam represented 
by the red bars, the same figure also shows how the sensitivity is increased at the DOF’s 
where the ABC is applied. Figure 26 represents the error localization at those places 
where the blue bars are raised. The stem plot with the red circle on the top indicates the 


magnitude and localization of the actual error. 


Figure 26 indicates, how the damage detection is correctly identified, showing 
how areas of strong sensitivity are very prone to detect the stiffness errors imposed on the 
beam. Additionally the increased sensitivity due to the pinned ABC at a DOF adjacent to 
the element where the stiffness error perturbation is applied showed the effectiveness in 


the damage detection. 
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Figure 25. Stiffness sensitivity distribution Case 1. 


MODE 1 


oO 2 4 6 6 10 12 14 16 18 20 
MODES 1:2 


Oo 2 4 6 8 10 12 14 16 18 20 
MODES 1:3 


“a 2 4 6 B 10 #12 #14 «#16 «618  ©620 
MODES 1:4 
“(= a6 5 a © oo S © «© «1 
a ai a ie a ie a a a 
0.1 
fn) 2 4 6 8 10 #12 14 #16 «#18 = 20 
MODES 1:5 





Figure 26. Stiffness error prediction Case 1. 
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2 Case 2. Stiffness Error at Element 5, and ABC Pinned at DOF 31 





Figure 27. Stiffness error at element 5, ABC pinned at DOF 31. 


Figure 28 shows the Stiffness sensitivity distribution along the beam represented 
by the red bars, this figure shows how the sensitivity is increased at the DOF where the 
ABC is imposed, which as for this case is located distant with respect to the element 


where the stiffness error is applied. 


Even though the sensitivity bars in the damaged element region are low, the actual 
error is located. This situation is due to the fact that the Base sensitivity distribution (no 
ABC’s applied) next to the fixed end of the beam is high, as it is shown in Figure 19, that 
is why the error prediction according to Figure 29 still preserves the damage when modes 


one through three up to modes one through five are retained. 


This case shows how even when the pinned ABC is applied at a DOF further 
apart from the stiffness damaged element, where no increase in the sensitivity distribution 
is achieved, it still can be able to retain the error perturbation, as long as the Base 


sensitivity distribution is high enough at that particular section of the model. 
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Figure 28. Stiffness sensitivity distribution Case 2. 
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Figure 29. Stiffness error prediction Case 2. 
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3. Case 3. Stiffness Error at Element 10, and ABC Pinned at DOF 19 





Figure 30. Stiffness error at element 10, ABC pinned at DOF 19. 


Case 3 is exactly another representation of case 1, where it can be seen in Figure 
31 the places of high sensitivity at the pinned DOF, which happens to coincide with the 


element among the actual stiffness error perturbation. 


This situation again predicts the error in a consistent way although not as accurate 
as in case 1. When modes one through three were retained, the error was picked up 
accurately, retained modes one through four and one through five bounded the error at 


element nine, which is still close to the place of the actual error location. 


One possible reason is due to the fact that areas of high strain energy 
concentration such as the fixed end of the cantilever beam are more likely to detect 
stiffness errors, thus the further away the damage element is from the fixed end of the 
cantilever beam model (area of high strain energy concentration), the more difficult the 
damage detection becomes. That is why the closest the damage element is to the fixed 
end, the better the damage detection will work due to the strain energy distribution. This 
is fortuitous, since the damage is more likely to develop near the fixed end, where the 


maximum stress concentration occurs (Panday, 1991). 
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Figure 31. Stiffness sensitivity distribution Case 3. 
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Figure 32. Stiffness error prediction Case 3. 








4. Case 4. Stiffness Error at Element 10, and ABC Pinned at DOF 31 





Figure 33. Stiffness error at element 10, ABC pinned at DOF 31. 


Case 4 is another example related to case 2, where Figure 34, shows how the 
sensitivity is increased next to the DOF where the ABC is applied, which as for this case 
is located distantly with respect to the element where the stiffness error is exerted, 
although, this same figure also shows how there is still some sensitivity exposure along 


the damaged element. 


The absence of blue bars in Figure 35 over the stem plot, shows how the error is 
never found, hence it confirms the notion that pinned ABC at a DOF further away from 
the element where the stiffness perturbation is applied is ineffective, in addition, since the 
Base Stiffness sensitivity or Base Strain Energy distributions are not really high at 


element 10, there are no possibilities that the error can be perceived. 


This case clearly demonstrates how in distant areas (DOF’s) where ABC is 
applied, the sensitivity is not very likely to be increased; hence no stiffness error is 


possible to be located. 
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Figure 34. Stiffness sensitivity distribution Case 4. 








MODE 1 


Oo 2 4 6 8 10 12 14 16 18 20 
MODES 1:2 


0.1 
a a ee ee: a a: ee 
0.1 


Oo 2 4 6 8 10 12 14 16 18 20 
MODES 1:3 


“oO 2 4 6 8 10 12 #14 «#16 «618 © «620 
MODES 1:4 
0.1 
rr se 
7 a a a a aE i a a ae 
0 2 4 6 8 10 #12 14 #16 18 = 20 
MODES 1:5 





Figure 35. Stiffness error prediction Case 4. 








5. Case 5. Stiffness Error at Element 15, and ABC Pinned at DOF 31 





Figure 36. Stiffness error at element 15, ABC pinned at DOF 31. 


Case 5 is analogous to cases | and 3 and concludes the set of representations of 
converged pinned ABC at the DOF location of the stiffness damaged element. This case 
again predicts the error in a consistent way, where the error detection is picked when the 
second mode is retained all the way up to retained mode five as it is depicted in Figure 


38. 


Cases 1, 3, and 5 prove the evidence found in Chapters IV and V where in 
particular for stiffness errors, when an ABC is particularly applied adjacent to the 
damaged element, the perturbation can be very well detected and retained, due to the 


increased sensitivity effect over the pinned ABC region. 
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Figure 37. Stiffness sensitivity distribution Case 5. 
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Figure 38. Stiffness error prediction Case 5. 


58 








6. Case 6. Stiffness Error at Element 15, and ABC Pinned at DOF 11 





Figure 39. Stiffness error at element 15, ABC pinned at DOF 11. 


Case 6 is the preceding and last description of cases 2 and 4, where the pinned 
ABC is applied at a DOF further away from the location of the stiffness damaged 


element. 


Figure 40 shows how the sensitivity varies arbitrarily with respect to the damage 
element, hence is very unlikely to have a reliable prediction of the error, as it is the case 
shown in Figure 41, nevertheless the error happened to get retained at modes one to four 
and modes one to five. Even though this was the situation aroused for this specific case, 
there is no consistent statement that this pinned ABC set up is going to accurately predict 


the stiffness error on the cantilever beam model. 


In summary, Cases 2, 4 and 6 ratify that when ABC’s are applied at a DOF further 
away from the element in error, there is not advantage of the increased stiffness of the 
fixed condition effect, therefore the sensitivity is not considerably raised; thus no 


stiffness error is likely to be found. 
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Figure 40. Stiffness sensitivity distribution Case 6. 
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Figure 41. Stiffness error prediction Case 6. 
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C. MASS DAMAGE DETECTION WITH ONE ABC 


ik Case 1. Mass Error at Element 5, and ABC Pinned at DOF 11 






Figure 42. Mass error at element 5, ABC pinned at DOF 11. 


Figure 43 shows the Mass sensitivity distribution along the beam represented by 
the magenta bars, this figure also shows how the sensitivity is decreased at the DOF 
where the ABC is applied. Figure 44 represents the error localization at those places 
where the blue bars are raised. The stem plot with the green circle on the top indicates the 


magnitude and localization of the actual mass error. 


As for this case, Figure 44 indicates, how no damage detection is identified, since 
no blue bars are visible along the twenty element beam model, showing how areas of low 


sensitivity are not going to detect the mass errors imposed on the beam. 


Since the mass perturbation error is directly correlated to the Kinetic Energy 
distribution of the beam, the mass error is a function of the velocity (inertia) of the model. 
Consequently the possibility of detecting mass errors at DOF’s where the ABC coincides 
with the element with the mass error perturbation are minimal. This artificial 
manipulation implies that the corresponding DOF is pinned to ground, therefore no 


displacement is very likely to occur. 
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Figure 43. Mass sensitivity distribution Case 1. 
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Figure 44. Mass error prediction Case 1. 
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2. Case 2. Mass Error at Element 5, and ABC Pinned at DOF 19 





Figure 45. Mass error at element 5, ABC pinned at DOF 19. 


Figure 46 shows the mass sensitivity distribution along the beam represented by 
the magenta bars, in the same figure it is also seen how the sensitivity is pushed away 
from the DOF where the ABC is applied, and how the sensitivity is increased at those 


DOF where the maximum displacement is present. 


Since element five is comprised within the area of maximum displacement, 
therefore good sensitivity is present at those DOF. The blue bars on Figure 47 depict how 
the error is well retained in magnitude and location when the first three modes are 
retained. This scenario shows how the mass perturbation errors are a function of the 


maximum displacement along the beam. 


Given that the mass perturbation error is correlated to the Kinetic energy 
distribution along the beam, the areas where displacement is evident, will also represent 


areas of Kinetic Energy presence; therefore the error is very likely to be predicted. 
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Figure 46. Mass sensitivity distribution Case 2. 
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Figure 47. Mass error prediction Case 2. 
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3. Case 3. Mass Error at Element 10, and ABC Pinned at DOF 19 





Figure 48. Mass error at element 10, ABC pinned at DOF 19. 


Case 3 is analogous to case 1, where the pinned DOF happens to overlap with the 
element with the actual mass error perturbation. Figure 49 clearly shows how the 


sensitivity is diminished at the DOF where the ABC is applied. 


The absence of blue bars on Figure 50 evidently proves again how no damage 
detection is identified. This case once again shows how areas of low sensitivity are not 
going to detect the mass errors imposed on the beam. Again, the lack of displacement 
effect due to the ABC pinned at a DOF adjacent to the element mass perturbation 


demonstrates how the mass error prediction is not likely to be successful. 


65 








Mode 1 


0.5 
oO 
0 2 4 6 3 10 12 14 16 18 20 
Mode 2 
1 
0.5 
0 
Zz 4 5 3 10 12 14 16 18 20 
Mode 3 
1 
0 
0 2 4 6 8 10 12 14 16 18 20 
Mode 4 
1 
0 
0 2 4 6 3 10 12 14 16 18 20 
Mode 5 
1 
oO 
oO 2 4 6 8 10 12 14 16 18 20 





Figure 49. Mass sensitivity distribution Case 3. 
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Figure 50. Mass error Prediction Case 3. 
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4. Case 4. Mass Error at Element 10, and ABC Pinned at DOF 31 





Figure 51. Mass error at element 10, ABC pinned at DOF 31 and DOF 13. 


Case 4 happened to be a very sturdy scenario on finding the mass damage 
detection. It was assumed that the pinned ABC that would generate one of the highest 
mass sensitivities measurements was that at DOF 31. Figure 52 shows the mass 
sensitivity distribution along the beam model for the mentioned pinned DOF. But Figure 
53 showed how the damage detection was never found, hence forcing to run the scenario 


for a variety of pinned ABC conditions. 


After several runs pinned ABC at either DOF 13 or DOF 15 (dark red triangle in 
Figure 51), proved to predict the mass error perturbation exactly in the same approach as 
it is depicted on Figure 55, where the magnitude and location of the error is found when 


the first three modes are retained and hold the error up to the five retained modes. 


DOF 13 happens to stay literally to the same distance with respect to element 10, 
as DOF 31 is, but the error prediction pattern was totally different when analyzed, where 
DOF 13 proved to do an acceptable error prediction, but the pinned ABC at DOF 31 did 


not work as efficiently as was expected. 


Although the evidence found in Chapters IV and V still holds true, that is, areas of 
higher displacements are very likely to predict better the mass errors, it is cumbersome, 
how it will not work for every case where the higher displacements are evident, making 


this situation a problem of major concern where mass error detection is of interest. 
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Figure 53. Mass error prediction case 4 at DOF 31. 
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Figure 54. Mass sensitivity Case 4 at DOF 13. 
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Figure 55. Mass error prediction case 4 at DOF 13. 








5, Case 5. Mass Error at Element 15, and ABC Pinned at DOF 31 





Figure 56. Mass error at element 15, ABC pinned at DOF 31. 


Case 5 is the last representation of cases 1 and 3, where the pinned DOF lies 
adjacent to the element with the actual mass error perturbation. Figure 49 clearly shows 


how the sensitivity is diminished at the DOF where the ABC is applied. 


Figure 50 proves again how no damage detection is identified. This case 
consolidates the reflection of how areas of low sensitivity are not going to detect the mass 
errors imposed on the beam. Again, the lack of displacement effect due to the ABC 
pinned at a DOF adjacent to the element mass perturbation demonstrates how the mass 


error prediction is never going to be successful. 
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Figure 57. Mass sensitivity distribution Case 5. 
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Figure 58. Mass error prediction Case 5. 
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6. Case 6. Mass Error at Element 15, and ABC Pinned at DOF 23 





Figure 59. Mass error at element 15, ABC pinned at DOF 23. 


Case 6 represented a very similar situation to the one experienced in case 4. 
Initially DOF 19 was selected as the pinned ABC condition to analyze, as it was assumed 
to generate a good Mass sensitivity distribution, but when the experiment was evaluated, 
it did not find the error accurately. Hence as in case 4, this case was run for a variety of 


pinned ABC’s conditions. 


The analysis was organized in two parts. The first with the ABC pinned to the left 
hand side of the mass damaged element, giving the boundary condition effect of a 
cantilever beam model. After several runs pinned ABC at DOF 23, was the best ABC 
configuration able to find the mass perturbation error, as it is shown on Figure 61, where 
modes one to three and one to four, retained the error in magnitude and location, mode 
five, lost the actual error detection, but still found some mass perturbation at elements 


fourteen and sixteen of the experimental cantilever beam. 


The second part of the analysis was with the ABC set up to the right hand side of 
the mass damaged element, giving the boundary condition effect of a clamped- clamped 
beam. The best error prediction was found at pinned ABC at DOF 35 and 37, which 
found the error at mode five only, the rest of the DOF did not find the error at all, and no 


error prediction pattern proved to exist. 
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Figure 60. Mass sensitivity distribution Case 6. 
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Figure 61. Mass error prediction Case 6. 
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D. STIFFNESS DAMAGE DETECTION WITH MULTIPLE ABC’S 


With the usefulness of the ABC’s proven for a single pinned ABC and one single 
damaged element, and using as guidance the different patterns found according to the 
different configuration of the ABC’s; it is inferred that by the use of more than one ABC, 
the system is going to become less underdetermined, therefore it is going to be able to 


better predict the error perturbations imposed in the cantilever beam model. 


The main goal is to configure the minimum number of ABC’s such that high 
enough sensitivity levels are attained throughout the structure, hence being able to have 
reliable damage detection capability at all times and at any beam location, considering 
that in a real life structure the location of the damaged elements are very likely to be 
unknown. This goal is going to be explored by the use of three and five ABC’s 
configurations along the cantilever beam, for one, two, and four damaged elements in 


order to test the ability of predicting the different stiffness errors in a reliable form. 


The following Figures are intended to represent damage detection for a different 
number and locations of elements with an imposed stiffness error, by the use of multiple 
ABC’s along the modeled structure. The sensitivity distribution graphs for the following 
scenarios were omitted, since it is already understood that the sensitivity levels are going 


to get increased a those DOF’s where the pinned ABC’s are applied. 


The pinned ABC’s are represented by the red triangles, the elements with the 
stiffness perturbations are represented in brown color. The stem plots with the blue 


circles on the top indicate the magnitude and localization of the actual stiffness errors. 
iF One Stiffness Damaged Element with Three ABC’s 


Figure 62 shows the cantilever beam representation for 10 percent stiffness 


damage at element ten, with a three ABC’s set up applied at DOF’s 15, 29, and 41: 
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Figure 62. Stiffness error at element 10, ABC’s pinned at DOF’s 15, 29, 41. 


The error prediction for the above configuration is the following: 

















Figure 63. Stiffness error prediction at element 10 with 3 ABC’s. 


Figure 64 shows the cantilever beam representation for a 10 percent stiffness 


damage at element five, and with the same three ABC’s set up: 
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Figure 64. Stiffness error at element 5. ABC’s pinned at DOF’s 15, 29, 41. 


The error prediction for the above configuration is the following: 

















Figure 65. Stiffness error prediction at element 5 with 3 ABC’s. 


2s Two Stiffness Damaged Elements with Three ABC’s 


Figure 66 depicts the same three ABC configuration, but now for a 10 percent 


stiffness damaged at elements 5 and 18: 


76 





Figure 66. Stiffness error at elements 5, and 18. ABC’s pinned at DOF’s 15, 29, 41. 


The error prediction for the above configuration is the following: 

















Figure 67. Stiffness error prediction at elements 5, and 18. 3 ABC’s. 


3. Four Stiffness Damaged Elements with Three ABC’s 


Figure 68 depicts the standard three ABC’s set up chosen for these experiments, 


for a 10 percent stiffness damage at elements 5, 10, 15, and 20: 
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Figure 68. Stiffness error at elements 5, 10, 15, and 20. ABC’s pinned at DOF’s 15, 
29, 41. 


The error prediction for the above configuration is the following: 

















Figure 69. Stiffness error prediction at elements 5, 10, 15, and 20. 3 ABC’s. 


4. One Stiffness Damaged Element with Five ABC’s 


Figure 70 shows the cantilever beam representation for 10 percent stiffness 


damage at element ten, with a five ABC’s set up applied at DOF’s 7, 15, 23, 33 and 41: 
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Figure 70. Stiffness error at element 10. ABC’s pinned at DOF’s 7, 15, 23, 33, 41. 


The error prediction for the above configuration is the following: 

















Figure 71. Stiffness error prediction at element 10. 5 ABC’s. 


Figure 72 shows the cantilever beam representation for a 10 percent stiffness 


damage at element five, and with the same five ABC’s configuration: 
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Figure 72. Stiffness error at element 5. ABC’s pinned at DOF’s 7, 15, 23, 33, 41. 


The error prediction for the above configuration is the following: 

















Figure 73. Stiffness error prediction at element 5. 5 ABC’s. 


a. Two Stiffness Damaged Elements with Five ABC’s 


Figure 74 depicts the five ABC’s set up, for a 10 percent stiffness damaged at 


elements 5 and 18: 
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Figure 74. Stiffness error at elements 5 and 18. ABC’s pinned at DOF’s 7, 15, 23, 33, 
4l. 


The error prediction for the above configuration is the following: 

















Figure 75. Stiffness error prediction at elements 5 and 18. 5 ABC’s. 


6. Four Stiffness Damaged Elements with Five ABC’s 


Figure 76 depicts the five ABC’s set up, but now for a 10 percent stiffness 


damage at elements 5, 10, 15, and 20: 
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Figure 76. Stiffness error at elements 5, 10, 15, and 20. ABC’s pinned at DOF’s 7, 
15, 23, 33, 41. 


The error prediction for the above configuration is the following: 














Figure 77. Stiffness error prediction at elements 5, 10, 15, and 20. 5 ABC’s. 
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E. MASS DAMAGE DETECTION WITH MULTIPLE ABC’S 


The following six scenarios represent the damage detection capabilities of 
multiple ABC’s for one, two, and four mass damaged elements. In the following Figures 
the pinned ABC’s are represented by the red triangles, the elements with the mass 
perturbations are represented in magenta color. The stem plots with the yellow circles on 


the top indicate the magnitude and localization of the actual mass errors. 
1. One Mass Damaged Element with Three ABC’s 


Figure 78 shows the cantilever beam representation for 10 percent mass damage 


at element eight, with a three ABC’s configuration applied at DOF’s 15, 29, and 41: 
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Figure 78. Mass error at element 8, ABC’s pinned at DOF’s 15, 29, 41. 


The error prediction for the above configuration is the following: 























Figure 79. Mass error prediction at element 8. 3 ABC’s. 
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2. Two Mass Damaged Elements with Three ABC’s 


Figure 80 shows the cantilever beam representation for 10 percent mass damage 


at elements five and fifteen, with a three ABC’s set up applied at DOF’s 15, 29, and 41: 





Figure 80. Mass errors at elements 5 and 15, ABC’s pinned at DOF’s 15, 29, 41. 


The error prediction for the above configuration is the following: 

















Figure 81. Mass errors prediction at elements 5 and 15. 3 ABC’s. 
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Figure 82 shows the cantilever beam representation for 10 percent mass damage 


at elements ten and eighteen, with a three ABC’s set up applied at DOF’s 15, 29, and 41: 





Figure 82. Mass errors at elements 10 and 18, ABC’s pinned at DOF’s 15, 29, 41. 


The error prediction for the above configuration is the following: 














Figure 83. Mass errors prediction at elements 10 and 18. 3 ABC’s. 
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3. Four Mass Damaged Elements with Three ABC’s 


Figure 84 shows the cantilever beam representation for 10 percent mass damage 
at elements five, ten, fifteen, and twenty, with a three ABC’s configuration applied at 


DOF’s 15, 29, and 41: 





Figure 84. Mass errors at elements 5, 10, 15, and 20. ABC’s pinned at DOF’s 15, 29, 
4l. 


The error prediction for the above configuration is the following: 

















Figure 85. Mass errors prediction at elements 5, 10, 15, and 20. 3 ABC’s. 
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4. One Mass Damaged Element with Five ABC’s 


Figure 86 shows the cantilever beam representation for 10 percent mass damage 


at element eight, with a five ABC’s configuration applied at DOF’s 7, 15, 23, 33, and 41: 





Figure 86. Mass error at element 8. ABC’s pinned at DOF’s 7, 15, 23, 33, and 41. 


The error prediction for the above configuration is the following: 














Figure 87. Mass error prediction at element 8. 5 ABC’s. 
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5. Two Mass Damaged Elements with Five ABC’s 


Figure 88 shows the cantilever beam representation for 10 percent mass damage 


at elements seven and eleven, with a five ABC configuration applied at DOF’s 7, 15, 23, 
33, and 41: 





Figure 88. Mass error at elements 7 and 11 ABC’s pinned at DOF’s 7, 15, 23, 33, and 
4l. 


The error prediction for the above configuration is the following: 

















Figure 89. Mass error prediction at elements 7 and 11. 5 ABC’s. 
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Figure 90 shows the cantilever beam representation for 10 percent mass damage 
at elements five and nine, with a five ABC’s set up applied at DOF’s 7, 15, 23, 33, and 
41: 





Figure 90. Mass error at elements 5 and 9. ABC’s pinned at DOF’s 7, 15, 23, 33, and 
41. 


The error prediction for the above configuration is the following: 

















Figure 91. Mass error prediction at elements 5 and 9. 5 ABC’s. 
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6. Four Mass Damaged Elements with Five ABC’s 


Figure 92 shows the cantilever beam representation for 10 percent mass damage 
at elements three, seven, eleven, and seventeen, with a five ABC’s set up applied at 


DOF’s 7, 15, 23, 33, and 41: 





Figure 92. Mass error at elements 3, 7, 11, and 17. ABC’s pinned at DOF’s 7, 15, 23, 
33, and 41. 


The error prediction for the above configuration is the following: 

















Figure 93. Mass error prediction at elements 3, 7, 11, and 17. 5 ABC’s. 
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Figure 94 shows the cantilever beam representation for 10 percent mass damage 


at elements five, nine, fifteen and eighteen, with a five ABC’s set up applied at DOF’s 7, 
15, 23, 33, and 41: 





Figure 94. Mass error at elements5, 9, 15, and 18. ABC’s pinned at DOF’s 7, 15, 23, 
33, and 41. 


The error prediction for the above configuration is the following: 

















Figure 95. Mass error prediction at elements 5, 9, 15, and 18. 5 ABC’s. 
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F. DISCUSSION OF RESULTS 


The sensitivity of the natural frequencies of the structure to changes in stiffness 


and mass were calculated from the finite element analysis. 


The twelve cases analyzed for the changes in stiffness; six for one ABC and one 
damaged element, and six for multiple ABC’s and multiple damaged number of elements, 


demonstrated a very specific behavior. 


For the one ABC and one damaged element cases, the highest probability of 
detecting the damage proved to be when the stiffness error was the closest to the fixed 
end of the cantilever beam, which happens to be the highest exposed part to stiffness 
failure due to the increased stress concentration in that area. The further the stiffness error 
perturbation was translated to the free end, the less likely it was detected, hence the use of 
the pinned ABC at a DOF adjacent to the element in error, highly increased the 


possibilities of detecting the damage. 


When the pinned ABC was applied to a DOF further away with respect to the 
location of the damage element, the error was never found, unless it was located very 
close to the fixed end of the cantilever beam, where the Stiffness Base sensitivity assisted 


in the detection of it. 


On the other hand, the multiple ABC and multiple damaged number of elements 
cases, significantly improved the damage detection capacity. For the 3 ABC’s 
configuration, the damage was successfully found up to two damaged elements, when 
four damaged elements were introduced, the system accurately predicted three out of the 
four damaged elements, and the fourth element in error although it was not detected, was 
very well bounded, as it is depicted in Figure 69. The five ABC configuration 
successfully found the damage imposed over the elements in error. Further experiments 
(not depicted in the present thesis) with more than only four damaged elements 
successfully retained the errors imposed over the system structure; this explains how due 
to the stable behavior of the ABC’s when to stiffness errors is concerned and making the 
problem less underdetermined with the incorporation of more ABC’s, the error is very 
likely to be found. 
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The twelve cases analyzed for the changes in mass; six for one ABC and one 
damaged element, and six for multiple ABC’s and multiple damaged number of elements, 


did not prove the same consistencies as the ones found in the stiffness errors. 


Since the mass errors are related to the inertia of the system, the pinned ABC’s 
were set up in such an approach that were able to produce a considerable displacement at 
the location of the damaged element, nevertheless, while this reflection holds true, not all 
of the different pinned ABC combinations at those DOF’ that produced higher kinetic 


energy seemed to detect the mass error perturbation successfully. 


For the one ABC and one damaged element cases, when the pinned ABC was set 
up to emulate the beam model as if it were a clamped-clamped beam, the error was 
sparsely detected and with no pattern defined; but the best error predictions were found 
when the ABC was configured to emulate the beam model as a cantilever beam, exerting 


the maximum inertia (displacement) at the damaged element. 


On regard to the multiple ABC’s and multiple damaged number of elements 
cases, the pattern found for the one ABC and one mass damage cases, stills holds the 
same. It was discovered that even when the configuration of more ABC helped to reduce 
the undetermined nature of the finite element problem, therefore being able to predict a 
few more elements in error, no specific pattern was found. Nevertheless it was 
recognized that the mass errors are a very strong function of the kinetic energy of the 
system. Figures 81, 85, 87, 89, and 93, show how the errors located at low kinetic energy 
regions were not detected, whereas, Figures 83, 91, and 95 show how the errors were 
accurately predicted when a smart combination of ABC produced areas of high kinetic 


energy at the regions where the damaged elements were located. 


On the following chapter an Optimization approach is develop in order to find the 
mass error in a cantilever beam. As it will be seen, the MATLAB built in optimization 


function, proves to work as a reliable technique. 
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VIL. OPTIMIZATION APPROACH 


The purpose of numerical optimization is to aid in rationally searching for the 
design variables to satisfy constraints (Vanderplaats, 1984). That is, to match the 
dynamic response of a finite element model to that of the experimental response of the 


system of interest. 


An optimization program was developed to predict the mass error 


detection/prediction on the cantilever beam model. 
A. OPTIMIZATION THEORY 


The general problem statement for a non-linear constrained optimization problem 


To minimize f(x) Objective Function 
Subject to: 
g(x) <0 (7.1) 
where j =l,m Inequality Constraints 
h, =0 (7.2) 
where k =1,p Equality Constraints 
KS hone (7.3) 
where i=1,n Side Constraints 


The design variables vector is represented: 
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c=. (7.4) 


The constraint equations are used to bound possible solution combinations to 
meet the requirements of the given situation. Most optimization algorithms require that 
an initial set of design variables, x° be specified. Beginning from this starting point, the 


design is updated iteratively. This process can be written as: 
x =x" +87 (7.5) 


where qrepresents the iteration number and s is a vector search direction in the design 


space. The scalar quantity a defines the distance that we wish to move in direction s. 
The updated values of x are used to calculate the value of the objective function for each 
iteration. The search direction is chosen to decrease the value of the objective function 
while staying within the constraints of the system. The process continues until the 


objective function converges to a local minimum and the optimal solution is localized. 
1. Modal Assurance Criterion (MAC) 


The modal assurance criterion is a scalar constant relating the causal relationship 
between two modal vectors. The MAC is used as a means to compare analytically and 
experimentally obtained vibrational mode shapes. The MAC will have a value between 0 
and 1. A value of 0 indicates that the two modal vectors are not correlated while a value 


of 1 indicates the modal vectors are correlated. The method of calculating the MAC for 


th 


comparing the i” mode of the analytical system to the i” model of the experimental 


system is: 
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(7.6) 


B. PROBLEM FORMULATION 


The optimization problem was analyzed according to the cantilever beam 
displayed on Figure 96, with the same beam specifications described on Table 1 from 
Chapter IV. The beam consisted of 42 elements, each 1 inch in length and DOF 1 was at 
the fixed end of the cantilever beam; this corresponded to FE model element quantity and 
length. A Series 336 FLEXCEL ICP accelerometer (serial number 10860) was threaded 
into position at node 42 of the beam and wired into Channel 2 on a DACTRON Focus 
front end digital signal processor (DSP). An excitation was applied by a PCB Series 086 
B03 impact hammer (serial number 269), which was wired into Channel 1 on the DSP. 
The accelerometer and force hammer were calibrated using a pendulous test; the 


sensitivity adjustment was applied in the set-up of RT Pro Focus 5.57 software. 





Figure 96. Experimental beam set up. 
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The experiment was performed in two segments: The analytical and the 
experimental. In the analytical, the Build2beams.m MATLAB code was used to generate 
the original finite element beam model set of original eigenvalues and eigenvectors, then 
the damage (mass perturbation) used during the real experiment, was simulated in 
MATLAB, generating the modified set of eigenvalues and eigenvectors; subsequently 
these modal parameters were evaluated into the MATLAB optimization program. The 
mass error perturbation was imposed on element 35 of the cantilever beam, with a 10% 
mass excess with respect to the total beam mass. The finite element model was analyzed 


from modes one to eight. 


The second (experimental) segment consisted of the actual measured data from 
the laboratory. It was created after tapping the beam shown in figure 62, producing the 
corresponding Frequencies Response Forces (FRF) by means of the softwares: 
DACTRON RT Pro Focus 5.57, and ME’scopeVES; node 42 remained the reference as 
the load cell roved from one node to the next allowing for the measurement of the 


response at each node. 


The modal data extracted out of the ME’scopeVES software comprises the natural 
frequencies, damping ratios and magnitudes and phases of the different mode shapes of 


the described cantilever beam with the mass perturbation added. 


1. The RT Pro Focus 5.57 software “Real-time” was configured to measure 3200 
spectral lines, 8192 points, with a delta T of 166.7s over the frequency range 0-2400Hz. 
The frequency range of 0-2400 Hz was chosen because it covered the first 10 modes of 
system and signal resolution was sufficient for data acquisition. The excitation signal 


proved to be clean and thus no window was used for data measuring. 


2. Channel set-up 

Channel | (Excitation) Channel 2(Response) 
Max Volts (mV): 0.1 0.3 
Quantity: Force Accel. 
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EU: lbf gn 
mv/EU: 8.67 104.383 
Coupling: ICP AC 0.7 Hz ICP AC 0.7 Hz 


Sensitivity Adjustment: 0 0 


3. Trigger set-up 

Source: Analog input 

Run Mode: Manual Arm every frame 
Input: Channel | 

Slope: Bi-polar 

Level (%): 1, Level (V): 0 

Pre/Post Points (-/+): -10 


Pre/Post Time (-/+): -1.67uUs 


4. Average set-up: 

Type: Linear 

Domain: Frequency 

Frames: 3 (Each node was excited 3 times and an average taken and saved.) 
Accept/Reject: Manual Accept/Reject every frame. 

(The user rejects double taps, under powered or overloaded signals.) 

5. Modal Coordinate set-up: 

Auto increment: ON 

Rove: Excitation 


Point increment: 1 
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Export: UFF text format, frequency response. 


Once the experimental modal parameters were extracted, a Complex to Real 
Conversion was performed on the mode shapes. Complex ones, as derived from analysis 
of test data, and real ones, such as are traditionally produced from theoretical analyses in 
the absence of a detailed knowledge of the nature, extent and distribution of damping 
(Ewins, 2000). The corresponding information was exported to MATLAB as a “.txt file” 
such that the corresponding experimental eigenvalues and eigenvectors could be read in a 


matrix form. 


The modal to spatial transformation was obtained by the use of the “Simple 
Method” (Ewins, 2000); through this method the modulus and the phase of the modal 
eigenvectors were obtained from the MEscopesVES software. The modulus being 
interpreted as the square root of the sum of the squares of the real and imaginary parts of 
the modal analysis data, and the phase was assigned to each node either as 0 or 180 
degrees, which represented the positive or negative value of the displacement part of the 
experimental eigenvector matrix. This process was performed within the MEscopeVES 
software once the shapes table was generated, by selecting “Complexity plot” from the 
“Display” drop down menu, obtaining a ‘Starbust plot’ — Argand Diagram (Ewins, 2000), 
then the modulus and phase table was automatically obtained after clicking on the 
“Normalize Data” button, now being ready to import the table as the “txt file” into 
matlab. Figure 97 depicts the Starbust plot obtained from the 10% mass error exerted at 


element 35 of the cantilever beam. 
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Figure 97. Starbust Plot- Argand Diagram. 


C. OPTIMIZATION APPROACH VIA MATLAB 


The actual optimization program was run with the use of a built in MATLAB 
function named “fmincon” which permitted the finding of the error perturbation on the 
beam. “fmincon” aids in finding a constrained minimum of a given Objective Function 
of several design variables starting at an initial estimate (xo). The fmincon algorithm 
employs the SQP implementation, which consists of three stages. The first stage updates 
the Hessian matrix via the quasi- Newton method, where the Hessian matrix is calculated 
using the Broydon- Fletcher Goldgarb-Shanno (BFGS) method, which is the most used 
with inexact line searches. The second stage performs the Quadratic Programming 
Solution, used to calculate the parameter changes, by the calculation of a feasible point 
and then the generation of an iterative sequence of feasible points that converge to the 
solution. The last stage essentially performs the line search and generates the merit 
function. The “fmincon” MATLAB set up used was: x = _ fmincon 


(fun,xo,A,b,Aeq, beg, lb,ub) 
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The “fmincon” needs as input arguments the following: 


e Objective Function (fun) 


Starting Point (xo) 


e Design Variable (DV) 


Equality Constraints (Aeq, beq) 


Inequality Constraints (A,b and lb,ub) 


1. Objective Function (OF) 


In order to prove the certainty and performance of the optimization program, in 


the mass- error detection analysis, three types of OF’s were implemented: 


a. Eigenvector Objective Function 


Taken from the Modal Assurance Criterion concept, which quantitatively 


checks the correlation between the measured and the calculated mode shapes, its 


arguments lie between 0 and 1; those elements with a value close to 1 indicate that the 


two modes under comparison are quite close. Thus the MAC was subtracted from one to 


produce a result that would indicate a better performance when minimized. 











b. Eigenvalue Objective Function 


Natural Frequencies Ratio 


P Phe la 
J(@)= > a 


=] 
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(7.7) 


(7.8) 


Cc. Eigenvalue plus Eigenvector Objective Function Combination 


Eigenvalues and eigenvectors were selected as the main core of the OF’s 


since they represent the modal response to the natural properties of the system. 


Pp 


J(g)= > 


1=1 








Ay = Ay 
Ay 


+]1 








(7.9) 








2. Starting Point (xo) 


Since it was always desired to iterate over the complete set of elements, the 


starting point consisted of a [42x1] unit value vector. 
3. Design Variables (DV) 


Being the eigenvalues and eigenvectors part of the natural properties of the 
system and since these are determined by the mass and stiffness of the structure, it was 
decided to use the density (g) as the design variable (DV). It was not intended to modify 
the elasticity modulus (E) of the beam, such that it could be used in future experiments. 
The design variable dv_rho was used under the modelo.m MATLAB program to 
assemble the finite element model with the corresponding density perturbation 


introduced. 
4. Equality Constraints 


No equality constraints were used for this experiment, hence its corresponding 


brackets in the matlab code were left blank. 
a Inequality Constraints 


A beam weight inequality constraint was set up in order to bound the weight of 
each element by a 10%, as well as the total amount of mass that could be added to the 


entire beam by a 10%. 
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D. RESULTS AND DISCUSSION 


The three Objective Functions were minimized subject to the following 
constraints: 
e System totally unconstrained. 
e System with limited amount of mass added to any one element (10 
percent). 
e System with limited total amount of mass added to the entire beam (10 


percent). 


The blue bars of the following figures represent the normalized damage detection 
distribution along the cantilever beam model. The top graphs of each figure represent the 
analytical damage detection, while the lower graphs represent the experimental damage 
detection distributions from the laboratory. The stem plots with the yellow circles on the 


top indicate the magnitude and localization of the actual errors. 
1. Modal Assurance Criterion (MAC) Objective Function 


Figures 98, 99 and 100 represent the analytical versus experimental comparison of 
the MAC objective function when it was used for the system totally unconstrained, for a 
one element mass constrained system and for the total mass constrained system, 


respectively. 


The eigenvector objective function proved to work very consistently, since the 
maximum value in the table was detected at element 35 where the actual error was 
located, although there was some error influence at the neighboring elements (34 and 36). 
The best damage detection/ error prediction performance was found at the total mass 
constrained condition, as it is shown in Figure 100, where the error is literally minimized 
and bounded at element 35, with no big influence on the rest of the elements of the 


cantilever beam. 
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Figure 98. Analytical- Experimental Unconstrained MAC O.F. 
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Figure 99. Analytical- Experimental Element Mass MAC O.F. 
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Figure 100. Analytical- Experimental Total Mass MAC O.F. 


2 Combined MAC and Eigenvalue Objective Function 


Figures 101, 102, and 103 represent the analytical versus experimental 
comparison of the MAC- Eigenvalue objective function when it was used for the system 
totally unconstrained, for a one element mass constrained system and for the total mass 


constrained system, respectively. 


This combined objective function proved to work almost as consistent as the 
MAC objective function did, yet again, the more constrained the problem became the 
more isolated the error was bounded as it can be seen in Figure 103. Theoretically 
speaking the more refined the objective function becomes, the better results are going to 
be found, in this case this combined objective function did not perform better than the 
MAC objective function alone, the most likely reason is due to the fact that the 
Eigenvalue objective function (equation 7.8) is not accurately predicting the damage 
therefore is not positively contributing to the combined objective function so it can build 


up a better performance. 
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Figure 101. Analytical- Experimental Unconstrained MAC/Eigenvalue O.F. 
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Figure 102. Analytical- Experimental Element Mass MAC/Eigenvalue O.F. 
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Figure 103. Analytical- Experimental Total Mass MAC/Eigenvalue O.F. 














3. Eigenvalue Objective Function 


Figures 104, 105, and 106 represent the analytical versus experimental 
comparison of the Eigenvalue objective function when it was used for the system totally 
unconstrained, for a one element mass constrained system and for the total mass 


constrained system, respectively. 


The use of the natural frequencies as objective function did not provide a good 
error estimate and seemed to work as the worst objective function out of the three 
selected. The main problem was the exceeded number of iterations with the program 
routine; hence this caused the program to be stopped, before meeting the desired 
tolerance conditions. The best error prediction was contained within the seventy percent 


of the actual error, not a very reliable situation. 
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Figure 104. Analytical- Experimental Unconstrained Eigenvalue O.F. 
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Figure 105. Analytical- Experimental Element Mass Eigenvalue O.F. 
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Figure 106. Analytical- Experimental Total Mass Eigenvalue O.F. 














The objective function that proved to work the best was the MAC; the combined 
MAC- Eigenvalue objective function did improve the results from just the natural 
frequencies objective function, but did not get better results with respect to the MAC 


objective function alone. 


The MAC objective function verified to work more efficiently than the 
eigenvalues objective function, proving the simplicity and low cost of the natural 
frequencies ratio objective function, but at the same time it confirms the disadvantages in 
the damage detection limitation, because the natural frequencies can not provide the 
spatial information about structural damage. On the other hand, the mode shapes 
objective function is more sensitive to damage, due to the strain energy correlation at that 
location (Jaishi and Ren, 2006). Further more the MAC, proved its usefulness in the 


mode shapes correlation, minimizing the objective function in question. 


110 


Theoretically speaking the more refined objective function should give the best 
damage detection results, but as it was proved with the combined MAC- Eigenvalue 
objective function, the reality not always matches with the theory, because it is not only a 
matter of having a more complete/ refined objective function, but it is also required to 
assure that the chosen parameters assure reliable results by themselves in order to 
improve the results when combined with other parameters and not to weaken the actual 
performance of the more complex/ refined objective function. Even when we had the best 
computer program, if the objective function is not refined enough, the results will not be 


satisfactorily found. 
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VII. CONCLUSIONS AND RECOMMENDATIONS 


Artificial Boundary Conditions have been verified to work as a method to find 
additional natural frequencies by using only one experimental data base, as opposed to 
reconfiguring the complete experiment such that more than the baseline frequencies 
could be found. Furthermore, the main purpose of this thesis was to accomplish the 
damage detection/ error prediction on sensitivity based finite element models by the 


judicious use of the Artificial Boundary Conditions. 
A. CONCLUSIONS 


1. The use of Artificial Boundary Conditions via the reduced order model 
proved to be an effective method to solve for the spatially incomplete real 
world structures. Hence reducing the disadvantages of the 


underdetermined problems. 


2. The sensitivity based updating governing equation proved to be a very 
useful tool in localizing the error/ damage, when the Artificial Boundary 


Conditions frequencies were added to the baseline frequencies. 


3. The stiffness sensitivity distribution was directly correlated to the strain 
energy distribution of the system; consequently the stiffness sensitivity is a 


function of the flexural properties of the system structure. 


4. The mass sensitivity distribution was directly correlated to the kinetic 
energy distribution of the system; as a result the mass sensitivity is a 
function of the inertia of the system when an excitation frequency is 


imposed over the system structure. 


5. Damage/ errors were very likely to be detected at areas of high sensitivity 


distribution of the structure. 
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6. 


10. 


The damage detection/ error prediction was consistently improved when 
the ABC’s were selected in a manner that increased sensitivity in the 


desired region. 


Stiffness errors proved a very consistent behavior. The highest damage 
detection probability was when the stiffness was the closest to the fixed 
end of the experimental cantilever beam. The use of the pinned ABC at a 
DOF adjacent to the element error, consistently improved the damage 
detection possibilities. When the pinned ABC was applied to a DOF 
further away with respect to the location of the damaged element, the error 
was less likely to be detected. The use of multiple ABC’s considerably 
improved the error detection capability of the FEM problem. 


Mass errors did not provide consistent evidence of damage detection. The 
only conclusion found was that since the mass errors are related to the 
kinetic energy of the system, the pinned ABC needs to be set up such that 
it is able to produce a considerable displacement at the location of the 
damaged element. Pinned ABC’s applied to mass based sensitivity 
systems tends to diverge the sensitivity, hence any ABC applied at a DOF 
adjacent to the mass damaged element, is less likely to find the error 
perturbation. The use of multiple ABC’s did not show an improvement as 


far as the error detection capability is concerned. 


The optimization program considerably improved the mass_ error 
detection performance of the cantilever beam model. The more refined the 


Objective Function, the more accurate/ minimized results are produced. 


The modal to spatial transformation in the Optimization routine was 
accomplished by the use of the “Simple Method”, although this method 
proved to work, it is believed that some more accurate methods would 
give a more accurate data transfer/conversion, as far as the softwares 


interface between the MEscopeVES — MATLAB is concerned. 
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B. 


RECOMMENDATIONS 


i, 


Find a method of approach to reliably predict the damaged mass elements 


over the beam model. 


Find or develop some interface routine to successfully pass directly the 
data generated between the software used for the finite element model 
simulation, such MEscopeVES, and the software used for the optimization 


computer program (MATLAB). 


Try more than only one physical parameter of the structure to be used as 
design variables in the formulation of the Objective Function, to allow 
isolation of all modeling errors while not requiring excessive iterations, 


hence obtaining an improved optimization result. 
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A. 


APPENDIX 


ABCRUNTHRU_RLA 


Yo KAKKKKKK KAKA KKK KK KKK ABCrunTHRU rla. mm *#ERRERRR RRA ARK KR KEK KEE 


This program calculates the condition number of the following 
sensitivity matrices used to calculate the DV (error prediction). 
1) Base system only 5 modes (underdetermi ned) 

% 2) ABC system 10 modes 
3) Base system 5 modes + 5 modes from ABC system 
The last systemis calculated 3 times. Once for modes 1-5, 
another for modes 6-10, and again for modes 11-15 


% This programis called from Build2Beams. m 


% Written by Constance RS Fernandez, Spring 2004 
% Updated by John R Mentzer, Spring 2007 


%| NPUTS 

| Re ee 

% icnt_oset 

% T_sens_ tot 

% vect_lamtot 


% OUTPUTS 


% intervel 

% mode 

% start mode 

% ten 

% dv_cal ABC - matrix 

% dv_calABCten - matrix 

% dv_cal_BasePlus - matrix 
% cond_ABC - matrix 

% cond _ABCten - matrix 

% cond basePlus - matrix 


% ----INITIALI ZATI ON------- 
abc = 0; 

ten = 1; 

intervel = 1; 

int_abe = 1; 


for aa = l:icnt_oset +1 % number of conditions (base + ABC) 
ac =1; % reinitialize for each ABC system 


for mode = 1:3 %3 sets of 
% (modes 1:5, 6:10, 11:15) 
startmode = abc + ac; 


dv_a = [startmode: start mode+4]; % modes 
dv_dl = [ten:ten]; % mode 1 

dv_d2 = [ten+tl:ten+l]; % mode 2 

dv_d3 = [tent2:ten+2]; % mode 3 

dv_d4 = [ten+3:tent+3]; % mode 4 

dv_d5 = [tent4:tent+4]; % mode 5 


dv_c =[ten:tent9]; % the first 10 modes of each ABC system 
% (modes 1-10 only) 


y 
Se modes 
5 


if dv_a == [1:.25*size(T_sens_tot,2)]; % if sensitivit 
thas only 5 rows than the modes used are all 5, el 
Yused are first five (base) and a set of selected 
wof ABC system for a total of 10 modes 
dv_b = [dv_a]; 

else 
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modes per ABC system (10 element beam 


matri x 


modes 


dv_b = [1:.25*size(T_sens_tot,2), dv_al; 


end 
%---Base System onl y---(underdetermi ned) 


% save DV calculates of as matrix for plo tting 

dv_cal ABC(:,intervel) = T sens _tot(dv_a,:)\vect clam Cott ay. a); hold on 
% Condition number of Sensitivity matrix used in DV Cal 

cond ABC(intervel,1) = cond(T_sens_tot(dv_a,:)); 


% cond ABC(intervel ae 
cond((T_sens_tot(dv_a,:))*(T_sens_tot(dv_a,:))'); %Ttrans)(T) 


% ‘mode 1 
% dv_cal_ABC1l = T_sens_tot(dv_dl,:)\vect_lam_tot(dv_d1); 
% %mode 2~ 
% dv_cal_ABC2 = T_sens_tot(dv_d2,:)\vect_lamtot(dv_d2); 
% %mode 3 
% dv_cal_ABC3 = T_sens_tot(dv_d3,:)\vect_lam_tot(dv_d3); 
% %mode 4 
% dv_cal_ABC4 = T_sens_tot(dv_d4,:)\vect_lamtot(dv_d4); 
% mode 5 
% dv_cal_ABC5 = T_sens_tot(dv_d5,:)\vect_lam_tot(dv_d5); 


% ---Base 
dv_cal_ Bas 


a des of ABC system teee 

e 

cond _basePlus(i 
Plu 
dv_ 


mo 
$(: ,intervel) T_sens_tot(dv_b,:)\ vect_lamtot(dv_b); 
ntervel, 1) = cond(T_ sens “tot(dv. b,:)); 
(inte rvel, 1) = 

1: ))*(T_sens_tot(dv_b,:))')); %Ttrans)(T) 


% cond base 


u 

i 
i. s 
cond(((T_sens_tot( b 


intervel = intervel + 1; 
ac =a_c +5; % five modes used at a time 
end % "mode" | oop 


% ---ABC system only ---- 


dv_cal_ABCten(:,int_abc) = T_sens_tot(dv_c, patina _lam_tot(dv_c); 

cond ABCten(int_abc,1) = cond(T_séns_tot(dv_c, 

%cond ABCten(int_abc,1) = 
cond(((T_sens_tot(dv_c,:))*(T_sens_tot(dv_c,:))')); %(Ttrans)(T) 


int_abc = int_abe + 1; 

abc = abc + 39; % advances to next ABC system, must change number "19" 
ae peuneat the number of DOF in beam. This beam had 10 elements thus 
0. 1 


ten = ten + 39; % advances to next ABC system 
end % "aa" loop 


Oy KARR KKK KKK KKK KEK KK END ABCrunTHRU rla. m #¥*KRRR HR RRRKK KKK EKE 
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B. ADDLUMPMASS_RLA 


Oy KAKA KKK KKKKKK KKK KK KK AddLumpMass_rla.m KEKKKKKKAKK KKK RR KKKKEK 


% ‘This script constru 
% which is added to t 
% Mass added 

H 


% Written by Prof J. 


cts a vector of lumped masses 

he diagonal of the BeamX mass matrix. 
to [mx] in Assembl e2Beams. m 

. Gordis 

% Inputs needed 


% num_elements 
% mx 


% Out puts 


% mass diag 


% mass node 

% mass coord 

% mass _ DOF 

% mx updated 

disp(’ ')idisp(' |); 

disp(' KKK KKK KKK KKK KK RK KKK KKK KKK KKK KKK KK KKK KKK KKK KAR KK KK KKK EK! 
disp(' **** Lumped mass addition to beams ¥xkE!) 
disp(' **** Lumpmass DOFs defined for UNRESTRAINED beams ****' ) 
disp(' KKK KKK KKK KKK KK RRR KKK KKK RK KKK KKK KKK KKK KKK KKK KKK KKK KK KK EK! 
% initalize 

if exist('mass_diag') == 0; % define and apply lumped mass vector. 
add_mass = 'n'; 

add_mass = input(' Add lumped masses to BeamX ? (y/n) ','s'); 


% Initialize vector to add to [mx] diagonal. 


mass diag = zeros(2*(num_elements+1),1); 


while add mass == 'y'; 
mass_node = input(' Node number for lumped mass ? '); 
mass_coord = input(' Translation or Rotation for lumped mass ? (t/r) ','s'); 


if mass_coord == ''t'; % Translational |umped mass 


mass DOF = 2 * mass_node - 1; 
elseif mass coord == 'r"; % rotational lumped mass 
mass DOF = 2 * mass_node; 
end 
mass_diag(mass_ DOF) = input(' Enter value of mass/inertia (in "Ibf-sec*2/in" '); 
% puts lumped mass on correct DOF 
add mass = input(' Add another lumped mass ? (y/n) ','s'); 
& can continue adding mass until 'n' is entered 


end; % End while loop 
end; % End exist('mass_diag') 


mx = mx + diag(mass_ diag); % Add |umped masses to [mx]: 


Oy KAKKKKKKKKKKK KKK KK KK END ADDLUMPMASS rla. m *¥#RRRKRKKKRRRRAEKKKKKE 
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C. ASSEMBLESENS_RLA 


Oy KAKA KKKKKKKKK KKK KK KK Assembl eSens rla. m  *k#EERRREKK KKK RRR KEKE 


x, T_sens_tot and 
relative frequency 
d Beam X 


% This program assembles the total sensitivity matri 
% total lam vector, vect_lamtot and assembles the 
% error between the natural frequencies of Beam A a 
% Written by Constance RS Fernandez, Spring 2004 


% Updated by John R Mentzer, Spring 2007 


n 


% | NPUTS 


% vect_lamx_oset 
% vect_lam_oset 
% vect_lam 

% T_sens_oset 

% T_sens 


% OUTPUTS 


% vect_ OSET 

% vect_lamtot 
% T_sens tot 
Se 

h freq X 

% rel_freqERROR 


vect _OSET = vect_lamx_oset - vect_lamoset; 

% |amx from actual beam with error oset, lam from FE model oset 

% Creating a vector of lamdifferences calculated (Lx-La) 

if vect OSET == 0; 
% when vector is empty at first, the total vector is equal to the 
%lamvector of BeamA, i.e., the first 19 values of vect _lam_tot are 
% the natural freq squared (rad*2/sec*2) of BeamA 


vect_lamtot = vect_lam 


else 
vect_lamtot = cat(1, vect_lam, vect_0OSET) 

end 

if T_sens_oset == 0; 

T_sens_tot = T_sens; 
else 
F T_sens_tot = cat(1, T_sens, T_sens_oset); 
en 


freq _OSET = sqrt(abs(vect_OSET))/2/pi 

% Natural frequency vector of Beam A in Hz 

freq OSETx = sqrt(abs(vect_lamx_oset))/2/pi; 

% Natural frequency vector of Beam X in Hz 

rel _freqERROR = freq_OSET. ie OSETx*100; 

% Relative error between Beam A OSET natural freq and Beam X OSET natural Freq. 


OAK ERR KKK KE REE ERE E END Assembl eSens rla. m #*RRERRRRRR KKK KERR KKK 
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D. 


BEAMPROPERTIES_RLA 


Oy KAKKKKKKKKKKK KK KK KKK BeamProperties Cla, M RRR RAKAKKRR RR RK KKK KES 


is the "props file" to load nominal beam data. 
programis called by BeamA_Prompt_crs to provide beam properties 
der to build Beam A. 


% This program was written by Constance Fernandez, Spring 2004 
% This progam modified by John R Mentzer, Spring 2007 

% Outputs 

Opa seh ece 

% depth 

% width 

% E 

% rho 


% total length 

% num elements 

% nominal El 

% nominal area 

% nominal density 


% Following are actual measurements from experimental set-up cantilever 


% beam 

depth = .504;% in z-dir (inches) 

ade a % in y-dir (inches) 
= e ' 
HE = 1.65e6; % | bf/sec%2-in (10e6- ksi ) 
%( 1bf/in*2 = 6894. 76Pa)-> E( I bf/in%*2) = ()Pa/6894. 76 

rho =0.110460934; %0.098; % | bf/in%3 
% T6 temper alloys require a 35-ksi tensile strength, 30-ksi yield 
% strength and a l0e6-ksi elastic modulus. Alloy 6061-T6 has 1.0 
% pct magnesium, 0.6 pct silicon, 0.3 pct copper and 0.2 pct chromium 
% lt has a 45-ksi tensile strength and 35-ksi yield strength.1 The 
% machinability of aluminum alloys are high (300) compared to titanium 
% (40). Aluminum alloys can easily be bent and provide easy loading and 
% unloading of parts. Also, aluminumis a highly conductive metal 


% compared to titanium 


% all measurement of distance are in inches 


total length = 42; 

num_el ements = 20 

nominal El = (width * depth*3 / 12) * E; 
nominal area = depth * width; % in%2 
nominal density = rho; % | bf/in%3 


Oy KAKA KKKKKKKKK KKK KK KK END BeamProperties Cla, M FRR RRREKKRRRR KEKE K KR ERS 
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E. BEAMSENSITIVITY_RLA 


Oy KAKKKKKKKKKKK KK KK KKK BeamSensitivity Cla. M  ERERRRRAAKKR RRR REA K KEKE 


i Revision history: 


% Ver. 1.0: 4/4/95 Basic frequency sensitivities 
% Updated: Spring 2004 Constance R S Fernandez 


OY RRR RRR RRR EERE RRR R RR ERR RRR ERR ERR R RRR R ERK ERE RRR RRR EERE REE K ERE 


% Program Description: 


% This program calculates mode frequency sensitivities as given by the 
% equation, 


% qw*2 11k] im 

% Be = {P}' * | > w2 ] * {P} 
% qDV qDV qDV 

% 

% where: w = natural frequency 

% P} = associated mode shape 


% DV design variable 
% The right side of the above equation is considered the addition of 
% the sensitivity matrix wrt El and the sensitivity matrix wrt mass 


% ‘The programcalculates the stiffness and mass matrix partials by finite 
% differences. That is, for example, the [k] matrix is assembled twice, 

% once in for the nominal beam parameters, and a second matrix is 

% assembled based on a small perturbation of element mass and/or El. 


% This program makes use of the beam data created by the program 
% ‘"“Build2Beams. m.' 

% 1) Resets BeamX mass and El data to be the same as BeamA data 
% 2) Enters a loop to create a sensitivity matrix (T) 


% In loop 

% a) <A small perturbation of mass is added to BeamX on ele.1, 

% b) The mass matrix is assembled for this mass-perturbed beam, 

: and the mass matrix partial derivative is calculated as 

0 

% q[ m] [m_perturb] - [m baseline] 

% ses = sss ss ee eee ee 

a qDV qDV 

9 

0 

% c) The first column of Sensitivity matrix is calculated using 

% 

% sens_mass = [phi_base]'*(-lam_base)*mdelta*[phi_base 

% 

A (Note: A column of T corresponds to the respective element on beam. ) 
9 

0 

% d) T loop starts again but with the small change on element 2. 

‘ e) Difference calcula ted and second column of T is calculated 

9 

0 

% 3) loop continues until all columns of T are calculated, corresponding 
i to a small change added to the respective element 

9 

0 

% 4) The procedure is identical for stiffness sensitivity. 

% Except: 

% sens El = [phi_base]'*k_delta*[phi_base] 

% 

% 5)Combine the two T's into one complete sensitivity matrix : 

% T_sens = [sens_mass,sens aie 

% In words, the first set of columns (equal to number of elements) of 
% the combined matrix is the sensitivity mass wrt mass changes and 
% the last half is wrt El changes. 


Oy RRR RRR RRR RRR R EER RRR RRR RRR RRR ER RE RR RRR RRR KERR RRR RRR ERE RRR REE K ERE ERE REE EEE 


% Inputs Needed 


Of, GS ed a.caace Beater 
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% mass_Ibls 

% El_Ibls 

% element_mass 

% element El 

% num_rbm 

% num_el ements 

% |amx (experimently measured nat freq of beam) 


% Programs called 


% Assemble2Beams crs; 

% BoundaryConditions_crs; 
% f modes 

% fOset_from Aset 


% Outputs: 


% num_modes 
% |am_base 
% phi_base 
% Sens_mass 
% sens El 
% T_sens 
% dv_cal 
% freq_base 
% vect_lam 


% Start code: 


format long; 


Oy KARR KKK KKK KKK RRR KKK KKK KKK RR KKK KKK RRR KKK KKK KR RR KKK KKK KKK KKK KKK KKK KK KK 
Oy KKKKKKKKRKKKK RRR ERRKKRKEK TNL TL ALI ZATION ¥¥ RRR E RRR K RRR RR RR KKK KKK KKK 
Oy RRR RRR RR EERE RR RR RRR ERR RRR R ERR RRR KERR ER RRR HERR ERR RRR R EERE EERE 


mass change 


1; % Percent mass change for sensit ne calculation. 
El_change Icula 


1; % Percent El change for sensitivity ca 


element_mass; % Copy properties to retain them. 


element_mass_orig 
element El; 


element _El_orig 


gaa for number of mode frequencies to generate sensitivities for: 
[ Seis e Meso oA Oe een ER eK AO OTe UN ear Aa oa amare ais A 


%num_modes = i : I ¢ 
num. modes = 39; % 39 is maximum number of modes in a 20 element cantilever 


% beam. This is hard coded for quicker run of simulation. 


input(' Enter number of modes for Sens calculations>> '); 


start_mode = numrbm + 1; % Skip the rigid body modes. 


Oy KKK RRR RRR KK KK RRR RR KK KKK RRR RRR KKK KR RRR KKK KK KR RRR KKK KK KKK KKK KK KKK KKK KK 


% KERR KEEKKKEKE MASS SENSITIVITY CALCULATION LOOP ***¥*## RRR KARR RHEEKE 
Oy KARR KKK KKK KKK RRR KKK KKK KKK RRR KKK KK RRR KKK KKK RRR RK KKK KK KKK KKK KKK KK KKK KK 


sens_mass = 0; % initialize Mass based sensitivity matrix 
if mass_Ibls ~= 0; % From Beam XPrompt as user inputs 
for icnt_dv = 1:numelements; % loop to create sensitivity matrix 


% Resetting BeamXx Peper eos to BeamA properties 
clement mass ,2) = element_mass(:,1); 


% each element, one at a time will have a change in mass 
element _mass(icnt_dv,2) = element_mass(icnt_dv,2) * .. 

(1 + mass_change/100); 
Assembl e2Beams crs; % Run script to assemble beams. 
BoundaryConditions_crs; % Apply boundary conditions. 
[lam_base, phi_base] = fModes(ka, ma); 


% ma is the basebeam without changes, 
% mx is beam with small change in mass added (mass change) 


% Form mass derivative matrices: 
mdelta = (mx - ma)/(mass_change/100); 
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% converts precent change into decimal amount 


% Mode freq sens loop: 
end_mode = start_mode + (num_modes - 1); 
row num = 0; % initialize loop 
for icnt_modes = start_mode:end_mode; 
row_num = row_num #41; 
sens _mass(row num,icnt_dv) = phi_base(:,icnt_modes)' *... 
-lam_base(icnt_modes) * mdelta ) *... 
phi_base(:,icnt_modes); 
% definition of mass sensitivity matrix 
end; % for "“icnt_modes" inner loop 


end; % End "for icnt_dv" outer loop for sensitivity calculations 
end; % End "if mass_Ibls ~= 0" 


Oy RRR RRR RRR RRR RRR R RR RRR RR EERE RRR RRR KERR RRR R ERK ERE RRR RRR KERR REE EEE 


Y KKKKKKKKKKAKKE EL SENSI TIVITY CALCULATION LOOP **#**¥¥# RRR KKKKKKK KKK 
Oy KARR KR KR KKK KR RRR RK KKK KR RRR RK KKK KR RRR KKK KK KR RRR KKK KK KKK KKK KKK KKK KK KK 


s El = 0; % initialize El based sensitivity matrix 
ET_Ibls ~= 0; % From Beam XPrompt as user inputs 


for icnt_dv = 1:numelements; % loop to create sensitivity matrix 


% ee BeamX properties to BeamA properties 
element _El(:,2) = element_El(:,1); 


% each e | 
element El 
(1 + ( 


have a change in El 
nt_dv, 2) Pee 


* 


i | 
cnt_dv,2) = element El (i 
“change/ 100) ); 


Assembl e2Beams_rla; % Run script to assemble beams. 


ement, one at a time w 
i c 
| 


| 
(i 
E 


BoundaryConditions rla; % Apply boundary conditions. 
[lam_base, phi_base] = f Modes(ka, ma); 

% ka is the basebeam without changes, 

% kx is beam with small sensitivity added 

matrices: 


ive 
El _change/ 100); 
ange to decimal value 


% Form El deriv 
k_ delta = (kx - ka 
% converts precent 


% Mode freq sens loop: 
end_mode = start_mode + (nummodes - 1); 
row num = 0; 


for icnt_modes = start_mode:end_mode; 
row Num = row_num #1; 

sens_El(row_num,icnt_dv) = 

k_ delta * phi _ base(:,i 

%definition of El sensitiv 

end; % for "icnt_modes" inner 


end; % End 
end; % if El_Ibls = 0; 


phi_base(:,icnt_modes)' *... 
cnt_modes); 

ity matrix 

loop 


"for icnt_dv" outer loop for sensitivity calculations 


El and mass properties back into arrays: 
element_El_orig; 
element mass orig; 


element El 
element mass 





% Copy element 


% cleans up workspace by clears all unimportant parameters 
clear element _El_orig element _mass_ orig end_mode start_mode 
clear row_numicnt_modes El_change k_delta mass_change mdelta 
clear ka kx ma mx Tcnt_dv 


% Assembles complete sensitivity matrix 
if sens_mass == 0& sens_El ~=0; 


T_sens = sens El; % resultant sensitivity matrix is equal 
% to El sensitivity matrix with only El changes if no inputs 
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% are given for mass changes 


elseif sens_mass ~= 0 & sens_El == 
T_sens = sens_mass; % resultant sensitivity matrix is equa 
% to Mass be aS matrix with only mass changes if no inputs 
% are given for El changes 


else 
T_sens = cat(2, sens_mass,sens_El);% else the resultant sensitivity 
fmatrix is teh combination of mass ey matrix in first set 
%of columns and El sensitivity matrix in the last set of columns 


end 


freqx = sqrt(lamx)/2/ pi 

freq base = sqrt(lam_base)/2/ pi 

% NOTE: % lamx = experiment measured natural freq of beam 
vect_lam = (lamx(1: num _modes)-lam_base(1:num_modes) ) 


OY RARER RRR RRR REE REE E END BeamSensitivity Pha. M FERRER RERR EERE KEE ES 
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F. 


BEAMSENSITIVITYOSET_RLA 


Oy KAKKKKKKKKKKK KKK KK KK BeamSensitivityOSET Pha. mM ERR RRKRK KERR RRR KKK KR KKK 


% Revision history: 
Of: ated ees we > ead eins od. o8. 0 


% Ver. 1.0: 4/4/95 Basic frequency sensitivities 
% Updated: Spring 2004 Constance R S$ Fernandez to create resulting 
% sensitivity matrix using lam vector of multiple ABC systems 


CH=}, =[0}»zRPPR Ee eee eee ee Ee 


% Program Description: 


% This program calculates mode frequency sensitivities as given by the 
% equation, 


% qw*2 {[k] 4[ m] 

% -- = {PH # [wees = WD =e] * {PH 
: qDV qDV qDV 

% 

% where: w = natural frequency 

% P} = associated mode shape 


% DV design variable 
% The right side of the above equation is considered the addition of 
% the sensitivity matrix wrt mass and/or El. 


% ‘The programcalculates the stiffness and mass matrix partials by finite 
% differences. That is, for example, the [k] matrix is assembled twice, 
% once in for the nominal beam parameters, and a second matrix is 

% assembled based on a small perturbation of element mass and/or El. 


% This program makes use of the beam data created by the program 
% "Build2Beams. m. " 

% 1) Resets BeamX mass and El data to be the same as BeaméA data. 
% %2) Enters a loop to create a sensitivity matrix (T) 


% In loop 

% a) A small perturbation of mass is added to BeamX on ele.1. 

% b) The mass matrix is assembled for this mass-perturbed beam, 

% and the matrix partial derivative is calculated as 

% 

i q[ m] [m_perturb] - [m baseline] 

(i ses SS ss ss ss ee we ee 

% qDV qDV 

% 

% c) The first column of Sensitivity matrix is calculated using: 
% 

f sens_mass = [phi_base -lam_base)*m_delta*[phi_base 

: [phi_b ]'*(-1 b )*m_delta*[phi_b ] 

0 

% (Note: A column of T corresponds to the respective element on beam.) 
% 

% d) T loop starts again but with the small change on element 2. 

% e) Difference calculated and second column of T is calculated. 

% 

% 3) loop continues until all columns of T are calculated, corresponding 
% to a small change added to the respective element. 

% 

% 4) The procedure is identical for stiffness sensitivity. 

% Except: 

; sens_El = [phi_base]'*k_delta*[phi_base] 

0 

% 5)Combine the two T's into one complete sensitivity matrix : 

% T_sens = [sens_mass, sens El. 

% In other words, the first set of columns (equal to number of elements) of 
% the combined matrix is the coo aa with respect to (wrt) mass 
% changes and the last half is wrt El changes. 


% El_Ibls 
% element mass 
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% element El 
% num_rbm 
% mx_beam 
% kx_beam 


% Programs called 


% Assemble2Beams crs 

% BoundaryConditions_crs; 
% f modes 

% fOset_from Aset 

% displacment Pl ot_OSET 


% Outputs: 


% num_modes 0 

% sens _mass0O, sens_El0 
% num_ 7 bmOSET 

% T sens _oset 

% dispX_fot, dispA_tot 
% accel -plot 
% oset Choice 

% accel ometer 

% BC, BCose, BCOSET 

% remaindof 

% ASETtot, OSET, OSETtot 
% inct_sens 


% mass change, El change 
% phi XPLOT, phi APLOT 

% maQ_base, kaQ_base 

% mxO, kx0 

% poe plotkx 

% |amaOSET, phiaOSET 

% |amxOSET, phizOSET 

% |lamxplot, phixplot 

% fad 

% mdelta0, k_delta0O 


% Start code: 


Oy KKK KKK KKK KKK RRR KKK KKK KKK RR KKK KKK RRR KKK KKK KR RRR KKK KK KKK KKK KKK KK KKK KK 
Oy KKKKKKKKRKKAKRRRRRKKKRKERK [NL TL ALI ZATION FR RRR E RRR R RRR RR RR KKK KKK KEKE 
Oy RRR RRR RRR RRR E REE RRR R RRR RRR EERE RRR RRR EERE ERR RERKE RE RRR RRR RRR ERE ER EE 


T_sens_oset = [] 
vect lam _oset 
dispX_tof =[ 
dispA tot = [| 
inct_sens = 0; 
icnt oset = 0; 


BCOSET = zeros(ndof, ndof); 
OSETtot = zeros(ndof, ndof); 
ASETtot = zeros(ndof, ndof); 
oset choice = 'n'; 


ndof % prints ndof to screen for user's reference 
accelometer = [3:2:ndof]; %This needs to change if you do not use a clamped free 


beam 

waccelometer = [3 11 13 15 17 19 21]; % use this line for convenience 
% in quicker calc n loops or use the next 2 lines. 

Maccelometer = in n what nodes are the accel. located','s'); %requests 
% user to input | ns of accelometers 

Yaccelometer = ev ',accelometer,']']); % converts string to vector of 


% labels 


—co°ow 


u 
p 
0 
a 


% saved for plotting accelometers in correct position. 
eae = floor(accelometer/ 2) +1; 
% graphical representation of accelometer locations 


oset_choice = input(' Do you want to use ASET and OSET (y/n)? ','s') 
% user can choose to use ASET and OSET calculations 
while oset_choice ~= 'n'; 
icnt_oset = icnt_oset a 
disp(' locations of Accel. ') % dl spl ays "location of Accel" on screen 
accelometer % displays the vector of location of Accel for user's 
Tees ene in choosing which DOF are to be pinned 
isp(' t 
disp(' Enter pinned DOF label(s)' 
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disp(' Use MATLAB vector format> 13 5:7 9 ') 

pinned = input(' >> ','s') 

pinned eval ([' le Pi ned,']']); % Converts string to vector of labels 
“differ 


BC(icnt_oset,:) pinned]; % Boundry conditions for cantilever beam 
% change line if d ent beamis used 


BCoset = fOset_from_Aset(ndof, BC(icnt_oset,:)); % gets OSET from ASET 


BCOSET(icnt_oset, l:length(BCoset)) = BCoset; % copies OSET into another 
% vector to be used 


% loop to find remaining ac ce One ners: 
for icnt = 1: length(pinn 
remain(icnt,:) = fin dt bl neeCLCHE} == accelometer); 


end 


remaindof = fOset_from Aset((length(accelometer)), remain) 
% remaining unpinned DOF 


% remaining accelometers 


aset = accelometer(remaindof); 

ASETtot(icnt_oset, l:length(aset)) = aset 

% omitted oset, i.e. pinned accelometers 

OSET = fOset_fromAset(ndof, aset); 

OSETtot(icnt oset, 1l:length(OSET)) = OSET; % copies OSET into another 
% vector to be used in====== 

inct_sens = inct_sens+l; 

mass change = 1; % Percentage mass change for sensit yety eae ' pal 
El_change = 1; % Percentage El change for sensitivity calcu 
element_mass_orig = element_mass; % Copy properties over to retain them 
element _El_orig = element El; % 


% Prompt for number of mode frequencies for which to generate sensitivities 
CTI eager erate Pace gio 0 ee ie OE eRe pe ope py ER aye pepe OT 
% use the following three lines for user's input for number of modal 

% frequencies for which to generate Ser tea: or use the fourth line 

% which uses the maximum number of modes for the defined system 


%disp ('max number of modes ') 

%l ength(oset) 

‘num _modesO = input(' Enter number of elastic modes for sensitivity 
calculations>> '); 

num_modesO = ndof - length(BC(i et, 


cnt ); max number of modes in system 
aeealaa numrbm + 1; % Skip i 


d body modes 


0 :) 
nik rig 
lizes two vectors to be used in ploting program 


zeros(ndof,num_modesO); %cantilever beam 
zeros(ndof, num. modes0): wcantilever beam 


se 
he 


Oy RRR RRR EERE REE RR RRR RR RRR ERE R RRR RR KER ERE RRR REE ERR EERE R ER EEE EEE 


% KKKKKKKKKKKKEKK MASS SENSITIVITY CALCULATION LOOP **#*¥*¥#¥ RR XXKKKKKK KKK 
Oy KKK RRR KKK KK RRR KKK KKK KR RRR KKK KK RRR RR KKK KK KR RRR KKK KKK RRR EKA KKK KKK KK KE 


sens_massO = 0; 
if mass_Ibls ~= 0; % from Beam XPrompt as user inputs 
for icnt_dv = l:num_elements; % loop to create sensitivity matrix 


% Resetting BeamXx ge ae to aaehe ee ees 
Sl enent mast 2) = element_mass(: 


% each element, one at a time will have a change in mass 
element _mass(icnt_dv,2) = element ii have 2 chang uaeeg 
(1 + mass _change/ 100); 


Assembl e2Beams_crs; % Run script to assemble beams. 
% Articial Boundry Conditions 

maQ_base = ma(BCoset, BCoset); 

% new mass matrix of the system defined by OSET 

mxO = mx(BCoset, BCoset); 


plotmx = mx_beam(BCoset, BCoset); 
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different 


1: pi nned(1, 2)- ); 
4,: 


1: pinned(1, 2)-3 
4,: 


)y 


); 


% resulting mass matrix with ABC used in plotting 
plotkx = kx_beam(BCoset, BCoset); 
% resulting- stiffness matrix with ABC used in plotting 


kaOQ_base = ka(BCoset, BCoset); 
kxO = kx(BCoset, BCoset); 


% lam (natural fr eee rad*2/sec*2), phi (mode shapes 

[| amaOSET, phi aOSET] f Modes(kaO_base, maO_base); 

% natural freq of new artifically bounded base system 

[| amxOSET, phi xOSET] f Modes(kx0, mxO); 

% natural freq/ modes of new artifically bounded system with 

% either El or mass changes added similiar to orginial calculations 
Ru eee ee = fModes(plotkx, pl ot mx 

Yresulting lam & phi of ABC system used in plotting 


% Mode shapes 
if Length( pinned) ==1 & pinned == 3 %the next DOF acts a little 


% because of the location on beamthus it was easier 
% to phone aw 7 pea 
hi 





phi APLOT(2 = phi aOSET(1:ndof-3, :); 
phi XPLOT( 2: ater 2, } = phixplot(1:ndof-3, :) 
elseif length(pinned)>1 & pinned(1, 1) ==3 
phi APLOT( 2: pinned(1,2)-2,:) = phiaOSET(1: pinned(1, 2)-3 
phi APLOT( pinned(1,2)-1:ndof-3,:) = phi aOSET(pinned(1,2 
elseif l|ength( pinned) <2 
phi APLOT(1: pinned-3, :) hi aOSET(1: pinned-3, :); 
phi APLOT( pinned-1: ndof - 2, as = phi aOSET( pinned-2:ndof-3,:); 
phiXPLOT( 1: pinned-3, : phixplot(1:pinned-3, :) 
phi XPLOT( pi nned- 1: ndof- 2:4 = phixplot(pinned-2:ndof-3,:); 
else 
phi APLOT(1: pinned(1,1)-3, :) = phia ager ey pinned(1, 1)- Pa 
ee COTE IMS ed(1,1)-1:pinned(1,2)-3,:) PAT AOSETUEL andl. 1 
Ge reimnrer eats te = phi aOSET( pinned(1, 2)-3: ndof- 
phiXPLOT(1: pinned(1,1)-3, :) = phixplot(1:pinned(1, 1)- 
EL cet OL nine es ,1)-1:pinned(1,2)-3,:) = phi baretearnaed: 1 
"phi XPLOT( pinned(1, 2)-1:ndof-3,:) = phixplot(pinned(1, 2)-3:ndof- 


end % 


end; 


num_rbmOSET = length(find(|lamaOSET < 1)) 
% find number of rigid bodies in new ABC 
start_mode = numrbmOSET + 1; % Skip t 


“system 
he rigid body modes. 
faO = sqrt(lamaOSET)/(2*pi); %natural freq of the ABC in Hz 


% Form mass derivative matrices: 

m.deltaO = (mxO - maO_base)/(mass_change/100); 

% converts percent change to decimal value 

% NOTE: mass_change needs to be the same change used in 
% orginial mass sensitivity matrix calculation 


% Mode freq sens | oop: 
end_mode = start_mode + (num modesO - 1); 
row_num = 0; 


for icnt_modes = start_mode:end_ mode; 
row_num = row_num #1; 
sens_massO(row num, icnt dv) = ph 
{-lamaOSET(i cnt _modes) * md 
phiaOSET(:,icnt_modes 
end; % end "for icnt _modes" inner loop 


T(:,icnt_modes)' *... 
Io Piste 


% End "for icnt_dv" outer loop for sensitivity calculations 
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feel 
)-2:ndof-4, 


displacmentPlot_0S 


% This program ass 
% sens _beam(dispX) 
if dispX_tot == 
dispX_tot = di 
dispA_tot = di 
else 
dispX_tot = ca 
dispA_ tot = ca 
end % end “if disp 


end; 


O KERR EEREKE ERE KEE 
% KKKKKKKKKKKKKK EL OF 
DE ae Naa 


sens E10 = 0; 
if ET_lbls ~= 0; 


for icnt_dv = I:n 


% Bee in 
element _El(: 


g 
,2 


me 
n 


e 
ic 
Ic 


Assembl e2Beams _ 


% Artifical Bo 
maO_base ma ( 
mxO = mx(BCose 


mx_be 
kx_be 


e = ka( 


S 
kx(BCose 


For two 


1: pinned(1,2)-3 
4,:); 


1: pinned(1, 2)-3 
he) 





i XPLOT( pi 


end 


% end "if mass_ 


APLOT( 1: 
APLOT( pi 


APLOT( pi 


i XPLOT(1: 
i XPLOT( pi 


ET % calls a program 


embles the beam displacement vector for 
and base beam(dispA) under ABC 


when the vector is empty 

% displacement vector used in plotting 
tot, displ); 

Geena 


Ibls = 0" 


KRKKKK KKK RRR RK KKK KKK KKK KK KKK KKK KKK KKK KKKKKKKK KKK 


NSIETIVIETY CALCULATION LOOP **##RRRRRKHREKRR KERR 
PEPE CCUCOOSCCCCCCOSO CC CCeCESCCCCeeeeLO CCC eee er. 


um_elements; % loop to create sensitivity matrix 
BeamX properties Ds BeamA properties 
) = element _El(: 


nt, one at a 
t_dv,2) = ele 
hange/ 100) ): 


time will have a rena: in El 
ment_El(icnt_dv,2) * . 


rla; % Run script to assemble beams. 


undary Conditions 
BCoset, BCoset); 
t, BCoset); 
am( BCoset, BCoset); 
am(BCoset, BCoset); 
BCoset, BCoset); 
t, BCoset); 
freq*2, rad*2/sec%*2), phi (mode shapes) 
OSET] = fModes(ka0_ base, ma0 base); 
OSET] = fModes(kx0,mx0): %séns info 
plot] = fModes(plotkx, pl ot mx); *plot info 
ed) ==1 & pinned == 3%Special case of DOF #3 pinned 
:ndof-2, :) = phi aOSET( 1: ndof- 3, :) 
ndof-2, :) = phixplot(1:ndof-3, :); 
inned)>1 & pinned(1, 1) ==3 
:pinned(1,2)-3,:) = phiaOSET(2:pinned(1, 2)-3,:); 
inned(1,2)-1:ndof-3,:) = phi aOSET(pinned(1, 2)-2:ndof 
4 inned) <2 %For Eee de ABC 
ae 35. 8 phiaOSET( 1: pinned-3, :); 
nied. ndof- 2,:) = phi aOSET( pinned-2:ndof-3,:) 
rateere 3y, 2) hixplot(1:pinned-3, :); 
inned-1: ndof- ae = phixplot(pinned-2:ndof-3,:); 
ABC's 
pinned(1,1)-3, :) = phia SUP ee: pinned(1, 1)- ye 
nned(1,1)-1: pinned(1,2)-3,:) PA aOSET{ pr aned( 1. 1)- 
nned(1,2)-1:ndof-3,:) = phi aOSET(pinned(1, 2)-2:ndof- 
pinned(1,1)-3, :) = phixplot(1: pinned(1,1)- 
nned(1,1)-1:pinned(1,2)-3,:) = phi barettarnaedils 1) - 
nned(1,2)-1:ndof-3,:) = phixplot(pinned(1, 2)-2:ndof- 
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num_rbmOSET = length(find(lamaOSET < 1)); 

start_mode = numrbm0SET + 1; % Skip the rigid body modes. 
faO0 = sqrt(lamaOSET)/(2*pi); %natural freq of the ABC, Hz 

% Form El derivative matrices: 

k_deltaO = (kxO - kaO_base)/(EIl_change/ 100); % in %/ 100 


% Mode freq sens | oop: 
end_mode = start_mode + (num_modesO - 1); 
row _num = 0; 
for icnt_modes = start _mode:end_mode; 
row_num = row_num #1; 
sens ElO(row_num,icnt_dv) = phiaOSET(:,icnt_modes)' *... 


k delta0* phiaOSET(:,icnt_modes); 
end; %end "for icnt_modes" 


end; % End "for icnt_dv" outer loop for sensitivity calculations 


displacmentPlot_OSET % calls a program 


% This poo ran assembles the beam displacement vector for 


% sens_beam(dispX) and base beam(dispA) under ABC 
if dispX_tot == 0; 
dispX_tot = displ; 
dispA_ tot = displa; 
else 
dispX_tot = cat(1,dispX_tot, displ); 
dispA_tot = cat(1, dispA_tot,displa); 
end % end “if dispX_tot == []" 
end; % end "if El_Ibls = []" 


% Copy element El and mass properties back into arrays: 
element El = element_El_orig; 
element mass = element maSs_orig; 


ement_mass_ orig end mode start_mode 


clear element _El_or | 
nt_ cnt_modes | bls num_ET_dv num_mass_dv 


g e 
clear row num icn vi 


i 
d 
% assemble of total sensitivity matrix for ABC System 
if sens_massO == 0 & sens_EIO ~=0; 

T_sensO = sens_El0; % no changes in mass 


elseif sens_massO ~= 0 & sens_El0O == 0; 
T_sensO = sens_mass0O; % no changes in El 


else 
T_sensO = cat(2, sens_mass0O,sens_E10); 
% changes in both mass and E 
end %end "if sens_massO == 0 & sens_EIO ~=0" 


% Builds complete sens matrix of all ABC systems 
if T_sens_oset == 0;% when matrix is originally empty 


T_sens_oset = T_sens0; 


else % after the matrix has some values 
T_sens_oset = cat(1,T_sens_oset, T_sens0); 


end %end "Tf T_sens_oset == []J" 
lamaOSET = |amaOSET(find(lamaOSET > 1)); %skips Rigid body modes 
vect_lamO = |amaOSET(1:num_modes0O); % nat freq of base beam under ABC 


builds a vector of natural freq of ABC systems 
if vect_lam_oset == 0;% when matrix is originally empty 


vect_lam_oset = vect_|am0; 
else% after the matrix has some values 
vect_lam_oset = cat(1, vect_lam_oset, vect_!am0); 
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end; 


aye aay vect_lam_oset == 0" 

isp(' ' 

oset_ choice = Hee Another cycle of ASET and OSET (y/n)? ','s'); 
% loop runs until "n" is inputed creating sens matrix & lam vector 
% for all ABC systems 

disp(' ') 


% end "while oset_choice ~='n'" 


Oy KK KKK KEK KKK KKK KEK KK BeamSensitivityOSET_rla.m KAKKKKKKKKKKKKKKEKEKR ARES 
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G. BOUNDARYCONDITIONS_RLA 


Oy KAKA KKKKKKKKKK KKK KKK BoundaryConditions rla.m KRKKKKKKK KKK KKK KKK KK KKK 


% Written by Prof Gordis 

% This script prompts the user boundary condition information 

% The script creates a vector of DOF (with ber pect to the unrestrained 

; saree and then extracts the rows and columns of the complementary 
0 . 


% ‘Script defines vector "free _dof_set" containing 
% (list of unrestrained dof 


% The boundary conditions are applied in this script. 


% Inputs needed: 


% free dof _set 
% updated ka, ma, kx, mx 
% icnt_dof 
% add dof 
% bc_node 
% bc_coord 
% bc_DOF 
% bc_boolean 
% all _dofs 
% restraint switch 
% KKK KKK KK ERK KK KK KARR KKK KK KR RRA KRK KAKA K KKK KR KKK KK KK ES 
% Start code 
if exist('free_dof_set') ==0; % Build free_dof_set vector 
disp(' Select a boundary condition set:') 
disp(' (1) Clamped-free') 
disp(' (2) Clamped- Cl amped' ) 
disp(' (3) Pinned-Pinned' ) 
disp(' (4) User-Defined') 
disp(' (5) Free-Free') 
BC_Choice = input(' >> Enter choice: ') 
if BC_Choice == 1; % Clamped-free =e 
free dof set = [3:ndof]; 
restraint switch = 'y' 
elseif BC_Choice == 2; % Clamped-Clamped =e 


free dof _set = [3:ndof-2]; 
restraint_switch = 'y'; 


elseif BC_Choice == 3; % Pinned- Pinned 


free dof set = [2:ndof-2 ndof]; 
restraint switch = 'y' 
elseif BC_Choice == 4; % User-Defined = = = se 
icnt_dof = 0; 
add dof = 'y'; 
whiTe add_dof == 'y'; 
bc_node = input(' Node number for restraint ? "0" to end: '); 
if bc_node == 
break 
end; 
be_coord = input(' Translation or Rotation ? (t/r) ','s'); 
icnt_dof = icnt_dof + 1; 
if be_coord == 't'; 
be_DOF(icnt_dof) = 2 * bc_node - 1; 
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end; 


ka 
ma 
kx 
mx 


elseif bc coord == 
bc_DOF(icnt_dof 
end;  % End if-els 


end; % End while add_dof 


bc_boolean = 
bc bool ean( bec 
alT_dofs = [1 
free _dof_set 

restraint swi 


elseif BC_Choice 


free dof set 
restraint _SWi 


end; 


ndo 


M peace 


== 5; % Free-free beam 


f, 


1); 


= 2 * bc_node; 
e block 


%®[111... icnt_dof] 


zeros(length(bc_DOF),1);% Put zeros in restrained dof 


% List of all dof 


SURE aR CRE POUCA ee Extract free dof 


[1: ndof]; 


tch ='n' 


% End if-elseif choice block 


% End exist block 


ka(free_ dof _set,free dof set); 


ma(free dof set, 
kx(free dof set 
mx(free dof set 


free dof 
free dof s 
free dof 


set); 
set); 
set); 


1 
1 
1 
1 


Oy KAKA KKKKKK KKK EEK END BoundaryCondi tions rla. mM ##ERRKKKKKKRKAKAKK KE KKE 
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H. BUILD2BEAMS_RLA 
% KEKKKKKKKKKKKKKK KEKE Bui l d2Beams rla.m KAKKEKKKKKKKKKKKKKKKK KK 
clear 
cle 
% Revision history: 
Of een eden eae oe 
% 
% Ver. 1.0: 9/22/94 Basic two beam assembly 
% 2.0: Added multi-element changes 
% 2.1 3/28/95 Added read/write to file, rebuild capability 
% 2.2 3/29/95 Added lumped mass additions 
% 3/10/04 Added Sensitivity matrices, error prediction, plots 


RAKKKKKKK KK KKK KK AKA KR KKK K KK KKK KKK KKK KEK KKK RRR KKK KKK KKK KKK KKK KKK KKK KKK 


Program Description: 


This program assembles the mass and stiffness matrices for 2 free-free 
beams, referred to as "BeamA" (analysis) and "BeamX" (experimental). The 
program can be run in several modes: 


"Build" mode: 


The user provides baseline data for BeamA, assumed to be a 
homogeneous, uniform beam Data provided: 


(1) Beam ra 

(2) Number of elements 

(3) Nominal El 

(4) Nominal cross-sectional area 
(5) Nominal weight density 


The program then prompts the user for instructions on how to modify 
"BeamA" data to arrive at "BeamX" data. The user can modify element 
masses, and/or element El values. The modification can be applied to 
either a single element, or range of elements, e.g. 

Modify single/range element mass values (y/n)? y 


If "y" is entered, the user enters the number of the element for mass 
adj ust ment: 


Enter element label(s) for mass modification: 1 
Use MATLAB vector format> 1 3 5:7 9 


Enter percentage mass change (+/- %) 
The user is prompted to modify another element or range of elements: 
Modify another element mass value (y/n)? y 


This process continues until the user enters an "n" for no change. 
This entire process can then be repeated for El adjustment. 


The program saves the beam definition data in a binary (.mat) file 
"beamdata" at the end of execution. 


The program can also be run in "Read" mode by entering an "r" at 
the initial prompt. 


Script Execution Path: 


Build2Beams_crs.m -- User executes this program. 
BeamA_Prompt_crs.m -- Prompts User for BeamA nominal beam data 
BeamX Prompt_crs.m -- Prompts User for BeamX modification beam data 
Assembl e2Beams_crs.m -- Called by Build2Beams, builds [ka] [ma] [kx] 
[mx], plots freqs. 
AddLumpmass_crs.m -- Prompts User for BeamX lumped mass addition 
BoundaryConditions_crs.m -- Prompts user for B.C.'s and applies them. 
PlotBeamModes_crs.m -- Calculate beam modes and plot frequencies 
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% BeamSensitivity_crs.m -- Calculate sensitivity matrix T-sens 

% BeamSensitivityOSET_crs.m-- Calculate sensitivity matrix using ABC 

% recorded Hcrs.m -- Calculates the nat. freq of BeamX with ABC 
applied 

% Assembl eSens_crs.m -- Assembles the sens matrices and calculates 
errors. 

% ABCrunTHRU. m -- Calculates the DV and cond number of matrix used 
% Saves data to "beamdata. mat" 


% Start code: 
Li) 


0 
Oy KKK KKK RK KKK KR RRR KK KK KR RRR RK KKK KR RR KKK KR KR RRR KKK KK RRR RRR KKK KKK RK RK KKK KK 


disp(' Building 2 beams from scratch...') 

BeamA Prompt _rla; % Prompt for BeamA Data: run prompt script 

BeamX Prompt _rla; % Prompt for BeamX Modification Data 

Assembl e2Beams_rla; % Run script to assemble mass and stiffness matrices 
la; % BeamX lumped mass vector construction and 


AddLumpmass_r 
% application 


kx_beam = kx; % saves the Beam X matrices without BC to be used later 

mx_beam = mx; 

BoundaryConditions rla; % Prompt for, and apply boundary conditions 

kx_beamBC = kx; % saves the Beam X matrices with BC to be used later 

mx beamBC = mx; 

be. Deenee = ka; % Saves the Beam A matrices with BC to be used later 
ma_beamBC = ma; 

PlotBeamModes _rla % Calculate beam modes and plot frequencies 

BeamSensitivity_rla; % Calculate sensitivity matrix T-sens 


BeamSensitivityOSET as Calculate sensitivity matrix using ABC 


recorded H rla; % Calulates the nat. freq of BeamX with ABC applied 
Assembl eSens ra; % Assembles the sens matrices and calculates errors. 
ABCrunTHRU_rTa; % Calculates the DV and cond number of matrix used 
FOM_rla; %Calculates the Figure of Merit for each prediction 
% Save Defining Parameters for Beams and plots 

disp(' ...saving beam data to file') 


save beamdata. mat 


Oy KAKA KKKKKKKKK KK KKK KK END Build2Beams rla. m #4 RRRRKEK KKK RRR KKK KK 
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I. 


DAMAGE_ PREDICTOR 


% % Mode error predictor 


% 
% 
% 
format long 
MLIPLL PTET LP ORGINAL METHOD/////7/007070701071T1 
calculates dv for base system summed 
for q=1:20 
4 dv_calculated_ base(:,q) = T_sens_tot(1:q,:)\vect_lam_tot(1:q); 
en 


%calculates dv for ABC1 system summed 
for t=40:59 

’ dv_calculated basel(:,(t-39)) = T_sens_tot(40:t,:)\vect_lam_tot(40:t); 
en 


calculates dv for ABC2 system summed 
for t=60:79 

; dv_calculated base2(:,(t-59)) = T_sens_tot(60:t,:)\vect_lam_tot(60:t); 
en 


calculates dv for ABC3 system summed 
for t=80:99 

4 dv_calculated base3(:,(t-79)) = T_sens_tot(80:t,:)\vect_lam_tot(80:t); 
en 


for t=100:119 


dv_calculated base3(:,(t-99)) = T_sens_tot(100:t,:)\vect_lamtot(100:t); 


end 


%For ABC System 1 and 2 166% %%%%%%% 109% %%%%% 100% WWD OY O 10 comentee 


Y% %%%%%%%%%%%%% 2 ABC System  — %%%%%%%%%%%%%%%% 0% %% %%%0 % 9% %0% 0% 10% %%%0% 
% T combol = zeros(9, 20); 


% T7combol(1:3,:)= T sens _tot(1:3,:); 

% T_combol(4:6,:)= T_sens_tot(40:42,:); 

4 T_combol(7:9,:)= T_sens_tot(78:80,:); 

0 

% vect_combol = zeros(9,1); 

% vect_combol(1:3,1) = vect_lam_tot(1:3,1); 
% vect_combol(4:6,1) = vect_lam_tot(40: 42,1) 
‘ vect_combol(7:9,1) = vect_lam_tot(78: 80, 1) 
0 

% for e = 1:9 


% a = T_combol(l:e,:)\vect_combol(1:e, 1) 
en 


% bbbb = (pcombol) 


% %%%%%%%%%%%% 3 ABC SYSTEM %%%%%%%%%%%%%%%%%% 10% 1% %% 
% a = reros( 12, 20); 

% T_combol(1:3,:)= T_sens_tot(1:3,:); 

% T_combol(4:6,:)= T_sens_tot(40:42,:); 

% T_combol(7:9,:)= T_sens_tot(78:80,:); 

: T_combol(10:12,:)= T_sens_tot(117:119,:); 

% 

% vect_combol = zeros(12,1) 

% vect_combol(1:3,1) = vect_lam_tot(1:3,1); 

% vect_combol(4:6,1) = vect_lam_tot(40: 42,1) 

% vect_combol(7:9,1) = vect_lam_tot(78: 80,1) 

% vect_combol(10:12,1) = vect_lam_tot(117:119, 1) 


% for e = 1:12 
% pcombol(e,:) 
% end 


% bbbb = (pcombol) 


= T_combol(l:e,:)\vect_combol(1:e, 1) 


WDWWWWWWWH% 5 ABC SyStSM %7%%%%%%%%%%%%0%%%%0%%%%0%%0%%0%% 
T_combol = zeros(18, 20); 

T_combol(1:3,:)= T_sens_tot(1:3,:); 

T_combol(4:6,:)= T_sens_tot(40:42,:); 

Le a ere T sens _tot(78:80,:) 

T_combol(10:12,:)= T_sens_tot(117:119,:); 
T_combol(13:15,:)= T_sens_tot(156:158,:); 


T_combol(16:18,:)= T_sens_tot(195:197,:); 


vect_combol = zeros(18,1); 

vect_combol(1:3,1) = vect_lamtot(1:3,1); 

vect combol(4:6,1) = vect_lam_tot(40: 42,1); 
vect_combol(7:9,1) = vect_lamtot(78: 80,1); 
vect_combol(10:12,1) = vect_lam_tot(117:119,1); 
vect_combol(13:15,1) = vect_lam_tot(156:158, 1); 
vect_combol(16:18,1) = vect_lam_tot(195:197,1); 
for e = 1:18 

HC ombo bea :) = T_combol(1:e,:)\vect_combol(1:e, 1) 
en 


bbbb = (pcombol)' 


% %For base and ABC System 16.46 %%%%%% 1% 9% %%%%%%%%%%%WWHOXOYO 10 comentee 


% %ABC 1 system 

% T_combol = zeros(6, 20); 

% T_combol(1:3,:)= T_sens_tot(1:3,:); 
Seep s: 6,:)= T_sens_tot(40:42,:); 

0 

% vect_combol = zeros(6,1); 

% vect_combol(1:3,1) = vect_lam_tot(1:3,1); 
4 vect_combol(4:6,1) = vect_lam_tot(40: 42, 1) 
0 

% for e = 1:6 

% pcombol(e,:) = T_combol(1l:e,:)\vect_combol(1:e, 1) 
% end 

% 

% bbbb = (pcombol)' 

% 

% %ABC 2 system 

% T_combo2 = zeros(6, 20); 

% T_combo2(1:3,:)= T_sens_tot(1:3,:); 

ae combo2(4:6,:)= T_sens_tot(78:80,:); 

0 

% vect_combo2 = zeros(6,1); 

% vect_combo2(1:3,1) = vect_lam_tot(1:3,1); 
paveck, comet’: 6,1) = vect_lam_tot(78: 80,1) 
0 

% for e = 1:6 


% pcombo2(e,:) = T_combo2(1l:e,:)\vect_combo2(1:e, 1) 
% end 

% ccccc = cond(T_combo) 

% cccc = (pcombo2)' 


%%%%%%%%%%%%%%%% %%0 20% % 0% 020% %% WWW WOWUWWDWD%Y 0 10 Comentee 
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J. 


% 


% Written by Constance Fernandez opt g 2004 
% Updated by John Mentzer, Spring 2007 
% This program plots the mode shapes (phi, |am) phi vs nodal position 
% of beam when ABC are applied. 
% Inputs 
| re 
% plotkx 
% kaO0_base 
% num modesO 
% phi XPLOT 
% phi APLOT 
% pinned 
% num_el ements 
Beagle 
%displ, displa 
%j ji, 9 
% ypos 
| era 
% displ = zeros(ceil(.5*size(plotkx,1)),num_modes0O) 
i nitilize disp vector and provides the first zero of the vector 
displ = zeros(ceil((20/38)*size(plotkx,1)), num modes0); 
Mi nitilize disp vector and provides the first zero of the vector 
% displa = zeros(ceil(.5*size(kaO_base,1)), num modes0); 
i nitilize disp vector and provides the first zero of the vector. 
displa = reros( cell (( 20/38) *si ze(ka0. base, 1)), num_modes0); 
i nitilize ee vector and provides the first zero of the vector. 
for jj = 1:ceil((20/38)*size(plotkx,1)); 
displ(jjtl,:) = phiXPLOT(2*jj-1, 1: num _modes0); 
% every other phi to give displacement at sequential nodes 
F displa(jj+1,:) = phi APLOT(2*jJ-1,1:num_modes0O); 
en 
% This loop normalizes the modes shapes to the tip modal displacement. 
if pinned == 41 % pt pinned is a special case, no new calcula 
ESET ECR 
ispla(:,:) = a(:, 
elseif fee Mates & pi nned(1, 2) ==41 %tip pinned with 2 ABC's 
displ(:,:) = displ(:,:); 
displa(:,:) = displa(:,:); 
else 
for g = 1: num _modesO 
displ(:,g) = displ(:,g)/displ(num elements +l, g); 
ental 1g) = displa(:,g)/displa(num_el ements+1, g); 
en 
i end 
mn en 
wypos = [l:1:num_elements+1l]; % Location of nodes used in plotting 
ypos = [1:1: num el ements +1]; % Location of nodes used in plotting 
vu mass_Ibls ~= []; 
0 
% for kk = lisize(mass_Ibls,1); 
% ff =0; 
% for JJ = l:length(find(mass_I bl s(kk,:)>0)); 
% ff = fftl; 
% posmkk, 2*]J-1) = mass _Ibls(k nk 
% posm(kk, 2*))) = mass_Ibls(kk, fy 
% end 
% end 
% 
% if kk == 1 
% posM = posm; 
% 
: else 
0 
% for uu = 1:kk-1; 
% posM = cat(2, posm(uu,:), posm(uutl,:)); 
% end 
% 
% end 


DISPLACEMENTPLOT_OSET_RLA 


RAKRKKKKKKKKKKKKK AKERS displacement Pl ot OSET rla. mM BHERKRRRRRHKK KKK KEKE EE 








posM = sort(posM(find(posM>0))); 
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tions are needed 


RARER KEKE RRR ERE RE EE 


_lbls ~= [] 
for iS Spree euaneieer 
for || = Ll:length(find(El_| bl s(kk,:)>0)); 
ff = ffl; 
pose(kk, 2*))-1) = El _Ibls(kk, ff) 
pose(kk, 2*)J) = El_I bls(kk, ff) +1; 
end 
end 
if kk == 
posE= pose 
else 
for uu = 1:kk-1; 
, posE = cat(2, pose(uu,:), pose(uutl,:)) 
en 


en 
osE = sort(posE(find(posE>0))); 
= Be ones (Ze DOSES) 41 ent ose 
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END displacement Pl ot_OSET_rla.m 


KKKKKKRKKKKKKR KKK KKK KKK KE 


K. 


FBEAMKM_RLA 


OY KAKKKKKKKAK KKK KKK KKK 


fbeamkm_rla.m 


% function [kbeam, mbeam] =f beamkm(|, ei, m) 
% Provided by Prof Gordis 


0 


Function [ kbeam, mbeam] =f beamkm(|, ei, m) 


KKKKKKKKKKK KKK KK KARR KK 


% This function returns the stiffness and mass matrices for 
% a simple 2-node beam element 


% Note 
9 


% Reference: R.D. 


% Outputs 


0 
% kbeam, mbeam, i, j 


kbeam=zeros 
mbeam=zeros 
% 


1 


4); 
4); 

0 

* 


= 





PPW PWN WNP 
Hin tt ton on nou 


4 


14; 
kbeam(j,i)= 
mbeam(j,i)= 
end 
end 


% 
i/1%3)*kbeam 
mbeam=(m/ 420. 0) *mbeam; 


% end function beamkm 


OY KKK KKKKKAKK KKK KKK KK 


m= rho * area * length = 


Cook, 


beami,j); 
J) 


37 
a 
o 
r-3) 
—% 


total 


Concepts and Applications of F.E. 


' 


END fbeamkm_rla.m 
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element mass 


Analysis 


RRKKKKKKKKKKKKKR KKK KARE 


L. 


FMODES 


YW REKKKKKKKAKK KKK KK KKK f Modes. m FERRER KKAKKKKKKEK KK KKK KK 


% function [lam, phi] =fmodes(k,m,num_to_print); 
% % Provided by Prof Gordis 


» % This program prints to the screen natural modes of system (phi). 


% % 
% % Usage: [lam, phi] =f modes(k,m,numto_print) 


% % ‘This function can be used with 1 to 3 arguments, as follows: 


% % [| am, phi] =f modes(a) : Produces modes of [a] with no print of freqs in 


% % [lam, phi] =fmodes(a,i) : Produces modes of [a] with print of 


"i" freqs in 


% % [lam, phi] =fmodes(k,m) : Produces modes of [m\k] with no print of freqs in 


% % This function returns a vector containing eigenvalues (rad/sec) *2 


% % and a matrix containing the mass normalized mode shapes. 
% % The mode information is sorted by frequency in ascending order. 


% % lf num_to_print > 0; tabular listing of numto_print freqs in Hz is 


% % lf num.to_print <= 0, no print. 


% % fNormalize 
% % Outputs 
% % phi 


% % error 


% 
% if nargin == 1; 
% [A] w/ no print request for freqs in Hz. 


% elseif nargin == 2 & size(m 1) == 1; % 
request for freqs in Hz. 

% 

1 normalization 


5] = sort(abs(diag(d))); 
ices); 
i v(:,indices), ‘one'); 


% elseif nargin == 2 & size(m 1) > 1; 
wi no print request for freqs in Hz. 
% 


% 


w/ print 


[k],[ ml 


% elseif HL == 3 & size(k,1) > 1 & size(m 1) > 1; % [k],[m] w/ print 


request for 
% 


reqs in Hz. 


% mass normalization 


% [v,d] =ei g( mk); 

% [lam, index] =sort{abs( di ag (d))); 

i [phi] =fNormalize(v(:,index),'mass', m); 

0 

% else 

% 

% num_to_print = -1; 

: Dre ErToE from fModes.m: Check input arguments.') 
0 

% end 

% 

% if numto_print > length(k); 

% numto_print = length(k); 

: end 

0 

% if nargin <3 & eee 2)==0 & k(L:length(k)/2,1:length(k)/2) == 
zeros(length(k)/2,length(k)/2); % Have [A] matrix 
% e = 1; % Eigenvalues are wn 

% else 

% e = 0.5; % Eigenvalues are wn%2 

% end 


% 

% disp(' '),disp(' ') 

% di Spl! ~wwww www nnnne r) 

% disp('Freqs in Hz.:') 
h u i 
% disp((lam(1l:num_to_print). *e)/2/pi) 
% di Sp(' wrwwr nw wnnnne *) 

% 

i end 


0 
Y % RAK KKKKKKKK KEK KK KKKE END f Modes, m *¥RRRRREAKKKRRKKEKKKKKE 
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M. FNORMALIZE 


YH RAK KKKKKKAKK KK KKK KKK f Normal ize.m KEKKKARKKEKKKKRKKRRKKKK KEKE 
% function [phi] = fNormalize(phi, method, m); 
y 


0 
% A Usage: [phi] = fNormalize(phi, method, m); 
0 


% % phi: matrix whose columns are to be (independently) normalized. 
OR method: String variable. The following choices are available: 
0 10 

% % "mass! Mass normalization 

% % ‘inf! Infinity normalization 
% % ‘one’ First element = 1 

% % ‘length’ Length = 1 

% % 

% % m: matrix used for normalization, i.e., phi'*m*phi = eye 

% % 

UO IO NP en ope a OS, oh ee 

% % 

% switch method 

% 

: case 'mass' % Mass normalization 

0 

% % disp('mass normalization') 

: phi = phi * diag(sqrt(diag((phi' * m* phi). %*(-1)))); 
0 

% case ‘inf! % Infinity normalization 
% % disp('inf normalization') 

% for icnt_cols = 1:size(phi, 2); 

% phi(:,icnt_cols) = 
phi(:,icnt_cols)/norm(phi(:,icnt_cols),inf); 

% end 

% 

% case ‘one! % First element = 1 

% % disp('one normalization' 

si phi = phi * diag((phi(1,:).%(-1))'); 

0 

% case ‘'length' % Length = 1 

% % disp('length normalization') 

% for icnt_cols = 1l:size(phi, 2); 

% phi(:,icnt_cols) = 
phi(:,icnt_cols)./norm(phi(:,icnt_cols),'fro'); 

% end 

% 

% end 


YW RAR KKKKKKKKK KK KKK KKK END fNormalize. m #eERRRERR REAR KK RRR KK EE 
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N. 


FOSET_FROM_ASET 


OY KKK KKKKKAK KKK KKK KKK fOset from ASet. mM FERRER KKEKKKKRK KKK KKK 


pine [oset] = fOset_from Aset(ndof, aset) 
0 
% Usage: [oset] = fOset_from_Aset(ndof,aset); 


% This function determines the comp) emeEN eT subset "oset" 
% froma set [1:l:ndof] and the subset aset =[x xx...] 


of DOF. Set j 
(proper subse 
=n 


% ndof: Total number 
% aset: Retained DOF 
% oset: aset U oset 

t') 


% 
% Provided by Prof Gordis 


nset = [l:ndof]; 


beled "nset 
[1:1:ndof]) 


io 
=> 


for icnt = 1: leng ) 
r indices(icnt) nset == aset(icnt)); 
en 


bool = ones(size(nset)) 
bool(indices) = ae: 
> 


ze(indices)); 
oset = nset(find(bool ) 


) 


Oy FARRER KKK KKK RRR KKK EK fOset_from Aset.m RRKKKKKKKKKKKKK KEKE KEKE 
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O. 


FSPRINGMASS2 


OY KKK KKKKKAKK KK KKK KKK fSpringMass2.m KKK KKKKK KKK KKK KKK KKKEK 


Uno [k, ml) =f SpringMass2(springs, mass, BC); 


: ‘ Usage: function [k, ml =fSpringMass2(springs, mass, BC) 

0 70 

% % This function script assembles the stiffness [k] and mass 

% % [m] matrices for an assemblage of springs. 

% % 

% % 

% % A linear chain of springs and masses is assumed. 

% % The number of springs is defined by the length of the vector 
‘springs 

: : and their values by the elements of ‘springs. 

Oo 10 

% % The number of masses is defined by the length of the vector 'mass' 
% % and their values by the elements of 'mass' 

% % NOTE: The number of masses must equal to the final number of 
active 

% % DOF (i.e., after BC's applied) 

% % 

% % Boundary conditions are specified by the vector 'BC.' This vector 
5 contains the DOF numbers which are to be restrained 

0 7/0 

% % For example, to build the following system 

% % 

% % .01 02 2015 

% % [--///1/--[m)--///1-- “ml re -[m 

% % 5 6 

% % 

% springs = [1 1 ]; 

% mass =[.01 .01 ]; 

% BC = [1] 

% % 

% % 

DoW Tack 45 cp nate tee toe gsp teh ce orp as fre ss gin dopa h tenth bbe ante tte Wecgts sepeetel 

% % 

% % BEGIN SCRIPT 


% 
s 


if length(mass) == (length(springs) +l) - length( BC); 


k = zeros(length(springs) +1,length(springs) +1); 
m = zeros(length( mass) ); 
% assemble stiffness matrix: 


rows = [0 1] 
for ispring = 1: length(springs); 


rows = rows + 1; 
addthis = [springs(ispring) -springs(ispring);-springs(ispring) 
Pere Chena nay: 
k( rows, rows) = k(rows,rows) + addthis; 
end 
i f ~isempt y( BC) ; 
peels = fOset from_Aset(length(springs) +1, BC); 
k( keep, keep); | 
end 
% assemble mass matrix: 
= diag(mass); 
else 


ae ee in fSpringmass2. Check # masses, springs, and BC"s.') 
return 


end 
Oy KKKKKKKKKKKK KKK KKKKK END fSpringMass2.m KEKKKKKKKKKKKKKEKKEKEKKE 
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systems 


% This code solves for the Kinetic Energy function for structural 


KINETICENERGY 


format long 


cle 


las 


magnitudes 

ied them 

ives, 
(1/L*2*rotleft(i, 


rs ( 
ti pli 
ivat 


corresponding hermitian Shape function second der 


num elements 


i,k) +2/L*3*displeft( 
/kK)*L*24displeft(i,k) 





modes required 


Oi % Orie See ee WD 
*%*e HS to ~ a  _ 
MmN SS BR ee 
CNG MONK LN LN YD 
Ce et i 
~~ Nn ES 


Nsw toe Oa ee 


base matrix and then mul 





Nox tox xe 


out of the phi 
+1 
+ 


ee eS a 
He ae eee ee 


Mme ee Lm ee 





ting the or 
1 
| ) 
+ 
( 
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( 
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+ 
r 


iT) ee 
Wt — cow See Se ee Oa et ee 
se MRE EI Ow HS ma — nw a 
aS +t Om — Ve — MVS Ves Vie 


Ss generated to select the corresponding eigenvecto 
n 


length/num_el ements; 


lacements) 


area 
%El ement Lenght 


% rho 


L 





my eigenvector matrix a 42 X number of 


QqQoo.-ll.- war aamasvanx AMmam 





0.016067021; 


= 10e6; 
total 

n 
r 
| 
d 
f 
t 

. n 
(i 
d 
d 
t 
d 
d 
t 
L 
d 
t 
d 
t 


%Given data 


eo Oo 


%A 


ss 


.*j prime 
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quared/ 2) 


area*rho*kin_| ams 


(nominal _ 


T 


is to be 


not 


with one ABC applied 
The 42 (ndof) 
equation 


modes required 
with the corresponding Neral chan second 


part of the strain energ 


orm, 


base matrix and then multiplied them 
ral 


rresponding hermitian shape function second derivatives 


m elements 


1 


pinned + 1:42]; %This for 20 elements 


(: 


out of the phi 


1 
ts 
ET 


ector matrix a 42 X number of 
n 
S 


length/num_el ements; 


6067021; 
y calculated in a polynom 


derivative 


Kinetic Energy calculation of the base beam 
%t i mes i 


sitffnes/mass perturbations. 
%"i abcprime" corresponds to the vee 
a 


KINETICENERGYABC 
format long 


%al read 


cle 


Q. 
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quared/ 2) 


rgy equation: 


in ene 
area*rho*kin_| ams 


(nominal _ 


end 
% Assembling the strai 


end 
Tabc 


MULTIABCS_MODE_ERROR 


R. 


Boe base system and the identified ABC. 


programis to graph the error prediction wrt number of 
s is done for 


is 
Thi 


zeros(78, 20) 


1:78; 
norm T(t,1:20)= T_sens _tot(t, 


%This is purely a plotting code 
t 


%The purpose of thi 
Ymodes retained, 


norm T 


for 


:),int); 


:)/norm(T_sens_tot(t, 


hold on 


a Mm 


or Prediction for Base system’ ) 


w <—— --MM—_—o., 


— F< KE a SS 
Qo ~——wwwToL Su SE&66_ 
BAw. GOO. voeeaas, 
co co2o—— a—— 


end 


-1)/2),0,'9%', ((BC(:,3)-1)/2),0,' gh’, ((BC(:, 3)- 


:, 3) 


lot (((BC( 
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S. NORMRUNTHRU_CRS 


Oy KAKKKKKKKKKKKK KKK KKK NormRUNthru crs. m  ERRRERERKR RRR KKK KKK KKK S 


ogram finds the NORM of the columns of sensitivity matrix and the 
he rows of the inverse of the a ae matrix. This was used 
correlation of the good prediction to the ABC system used. It 


% also plots the information in helpful graphes 


s 

= 

oO 

 ~ had 

= 

°o 
Qa 


g 
t 
a 
t 


% This program was written for a system of 19 natural freq. Set of 5 modes 
% were used in each ABC system, i.e.,, modes 1-5, modes 6-10,or modes 11-15 
% This accounts for the 3 sets of modes per condition as listed belowin 

% the "for" |oop modeN = 1:3. This program compares the use of the first 5 
% modes of the base system and one set of five modes of the ABC to that of 
% of just 10 modes of the ABC. Notice that the sensitivity is a 10x10 
savers matrix indicating that only mass or El changes, not both were 

% made. 


% This programis not part of Build2Beams_crs.m program It is run 
% separately. 


% Written by Constance R S Fernandez, Spring 2004 


% cond _basePlus 

% icnt_oset 

% T_sens_tot 

% EL_Ibls, mass _Ibls 

% dv_mass, dv_ET, dv_cal_ABCten 


% Out puts 


% abcN, countN, a_cN 

% modeN, startmodeN, start modeNT 

% bb, t, tINV, cc, T, TINV, vv, tt 

% model abel NORM 

% norm_vectT, norm_vecTABC, norm_vecTinvABC 
% normc 

% baseN, baseABCN, abc_conN, abc_conTN 

% baseABCNten, abc_conNten, abc_ConTNten 


% ----Start program. - 
BASE = int2str(cond_basePlus(1)); % cond no. of the base line system 
BASET = sprintf('Base[1:5] Cond = %s', BASE); % used for plotting 
% ======/nitialization=====% 
abcN = 0; 
countN =0; 
% =======Calculations of NORM vectors======% 
for on ey icnt_oset +1 % number of conditions (base + ABC 
a_cN = 1; 


for modeN = 1:3 % 3 sets of modes per boundry condition 
startmodeN = abcN + a_cN; % the beginning mode number of each set 


% indicates the use of the first 5 modes of base system + 5 modes 
an of ABC pean 
= [1:5, startmodeN: start modeN+4] 


g of modes for plotting 
el NORM = int2str(a_cN: a_cN+4) 


a_cN = a_cN+t5; wadvances to the next set of modes 
lus 5 s of ABC 

bb,:); % 10x10 matrix 
ns_tot :)):% 10x10 matrix 


% base system 
t = T_sens Pate 
tINV = inv({T_se 


% first 10 modes of ABC solo 


startmodeNT = abcNtl; % the be i ung mode number of each set 
cc = [startmodeNT: start modeNT 
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T = T_sens tot(cc,:);% 10x10 matrix 
TINV = inv[{T_sens “tot(cc, :));% 10x10 matrix 


% for loop for NORM of columns and rows of inv(sens matrix) 
for vv = 1:10 % 10 rows, 10 columns 
% base + ABC system 
normvecT(vv,countN+tmodeN) = norm(t(:,vv)); % columns 
norm vecTinv( vv, countNtmodeN) = norm(tINV(vv,:)); % rows 
% for 10 modes ABC solo 
norm vecTABC( vv, countN+modeN) = norm(T(:,vv)); % columns 
norm_vecTinvABC(vv, countN+modeN) = norm(TINV(vv,:)); % rows 


end % vv loop 
end % ModeN loop 
abcN = abcNt+19; % advances to the next ABC system 
countN = countN+3; % counts up each set of ABC 
end % count loo 
normC=4; % initialize for plots 


0 % A as (30-40) plots 6 graphes per figure 


( 
m vect(: , nor mc) ) 
‘Norm col Tsens, ABC Modes [1:5]'); 


ot (3,2, 3) 
orm_vecTinv(:, norm) ) 


title ('Norm row TsensI NV, ABC Modes [1:5]'); 


subpl ot (3, 2, 2) 

bar( norm vecTABC(:, normC) ) 

title (‘Norm col Tsens, ABC Modes [1:10]') 
subplot (3, 2, 4) 

bar(norm vecTi nvABC(:, nor mC) } 

title (‘Norm row TsensINV, ABC Modes [1:10]') 





ah ot(3,2,5) 

% plotting error prediction 

% using [1:5] modes of ABC system + [1:5] modes of oe system; 
baseABCN = bar cdy. cal asePbuel cuernich 3, r');hold on 

% plotting error aco using [1:5] modes of Base system; 
baseN = bar(dv_cal_ABC(:,1),.25,'b'); 

abc_conN = int2str(cond_basePlus(normC)); % cond no. for legend 
abc_conTN = sprintf('Base[1:5]+ABC[1:5] Cond = %s', abc_conN); 


grid on 
legend([baseABCN, baseN], BASET, abc_conTN), hold on 


% ee actual error 


if El_Ilbls ~=[] & mass_lbls ==l 


=[] 

stem(mass_|bls, dv_mass, = ,'filled'); hold on; 
stem( mass | bls, dv_mass, ) 
Elplot = El_Ibls+10;hold on 
stem(El_Ibls, dv_El,'c','filled');hold on; 
stem El Ibls, dv_El,'k') 

elseif mass _Ibls ~=[] & El_Ibls ==[] 
stem(mass_| bls, realest a a a i on; 
stem(mass_I bls, dv_mass,'k') 

else 
stem El _Ibls, dv_El,'c','filled');hold on; 
stem(El I bls, dv_El,'k') 

end %if El_Ibls ~=[] & mass_Ibls ~=[] 





title (sprintf('Error, Base [1:5] + ABC [1:5], pinned NODE # %d', (tt+1))) 


Pa eee 
% plotting error prediction using [1:10] modes of ABC system; 
baseABCNten = bar(dv_cal_ABCten(:,normC));hold on 
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abc_conNten = int2st 


r(cond ABCt en( : 
abc_conTNten = sprintf('ABC[1:10] %s', abc_conNten); 
ee on 
egend([baseABCNten],abc_conTNten), hold on 


% plotting actual error 
if El_Ilbls ~=[] & mass_lbls =| 





] 

stem(mass_|bls, dv_mass, a) ,'filled'); hold on; 
stem( mass |bls, dv_mass, ) 
Elplot = El_Ibls+10; hold on 
stemEl _Ibls, dv_El,'c','filled');hold on; 
stem El _Ibls, dv_El,'k') 

elseif mass _Ibls ~=[] & El_Ibls ==[] 
stem(mass_| bls, ayes Seay d'); hold on 
stem(mass_|bls, dv_mass,'k') 

else 
stem El _Ibls, dv_El,'c','filled');hold on; 
stem(El _Ibls, dv_El,'k') 


end % El_Ibls ~=[] & mass_Ibls ~=[] 


title (sprintf('Error, ABC Modes [1:10], pinned NODE # %d', (tt+t1))) 


no 


rmC = normC+3; % advances to the next ABC system 
end % tt 


= 1:10 for plotting graphes 


Oy KKK KKK KKK KKK KEK KK END normRUNthru crs. m ¥¥RKRRRRRRRK KKK KKK KKK 
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FF PLOTBEAMMODES_CRS 


% KEKKKKKKKKKKKKKKKE KE Pl ot BeamModes crs.m KEKKKKKKKKKKKKKK KK KKKKE 
h Calculates natural frequencie 

% Calculat t | freq S 

% Provided by Prof Gordis 


% Inputs needed: 


% f Modes 
% Outputs: 
% lama, phia, |amx, phix (without rigid body modes) 


% num_rbm 
% phia_plot, phix_plot 


eee 
Calculating modes for each beam...plot frequency comparison') 


% Get modes of each beam: 


[| ama, phi a] =f Modes( ka, ma); 

[| amx, phi x] =f Modes(kx, mx); 

Yused to plot the mode shapes org BC before ABC 
phia_plot = phia; 

phix_plot = phix; 


% Set any rigid body mode freqs to zero: 


num_rbm = length(find(lama < 1)); 


sprintf('Number of Rigid Body Modes Found: %2i', numrbm) 


disp( ' Removing rigid body mode frequencies from vectors...') 
lama = |lama(find(lama > 1)); 
lamx = lamx(find(lamx > 1)); 


RARER KKK KERR ERE EE END PlotBeamModes_crs.m KKKKKKKKKK KKK KKK KKK KKK 
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U. PLOTTINGBARS_CRS 


RAKRKKKKKKKKKKKKKA KEE pl ottingBARS CES. Mm RERKKKK KKK KKK KK RK KKK K 


To be used with Build2Beams.m and 


This program plots 9 graphes per figure. The first columns of 3 graphes 
are the mode stapes of the ABC system used in error prediction. The next 
column of 3 utee es are the error prediction using only 5 modes of ABC 
system. The Tast column of 3 graphes are the error predictions using the 
first 5 modes of base system plus 5 modes of the ABC system The row 
represent modes 1-5, middle row: modes 6-10, last row: modes 11-15. Each 
of the error prediction graphes also have the base only prediction and 
the actual error plotted for easy reference 


Written by Constance R S Fernandez, Spring 2004 


cond basePlus 

FOM ABC5per, FOM PLUSper 
icnt_oset 

modeshape 

rel_freqERROR 


ypos 
El _Ibls, mass_Ibls 
dv_mass, dv_ET 


% Out puts 

| 

% BASE, BASET 

% FOMBASE, FOMABC, FOMPLUS 

% intervelp 

% ER, barp, shape, error, a_cp, modep, ap, 
% modelabelp, FOMABClabelp, FOMPLUSI abel p 
% abc_con, abc_conT 

% ABC, base 

% beeen plus_conT 

% base, plus, El_plot 

BASE = int2str(cond basePlus(1)); 

FOMBASE = int2str(FOM ABCS5per(1)); 

BASET = sprintf('Base Cond = %s, FOM = %s', BASE, FOMBASE) 
intervelp = 3; 

modeshape = 1; 

ER = 1; 

for barp = l:icnt_oset 


figure(barp+l0) % figures 11-20 
format bank 


%shape = [modeshape: modeshape+20]; %This number is equal to number of elements. 
shape = [modeshape: modeshape+20]; %This number is equal to number of elements 
error fia es freqERROR( ER: ER+15)*100)/ 100; 
a_cp = 
for modep = 1:3 %3 sets of modes per boundry condition 
[a_cp: a_cpt+4]; %modes 
PREL errorl = Iint2str(error(a_cp)); 
Errorlabelp = sprintf('Rel error = %s', REL_errorl) 
REL_error2 = int2str(error(a_cp+l)); 
Errorlabelp = sprintf('Rel error = %s', REL_error2) 
REL_error3 = int2str(error(a_cpt2)); 
Errorlabelp = sprintf('Rel error = %s', REL_error3) 
REL_error4 = int2str(error(a_cpt3)); 
Errorlabelp = sprintf('Rel e@rror = %s', REL_error4) 
REL_error5 = int2str(error(a_cp+4)); 
Errorlabelp = sprintf('Rel error = %s', REL_error5) 
modelabelp = int2str(a_cp:a_cpt4) 
FOMABC = | nt2str(FOM ABC5per(intervel pt+modep) ) 
FOMABCI abel p = sprintf('System FOM = %s FOMABC) 
FOMPLUS = int2str(FOM PLUSper(intervel pt+modep) ) 
FOMPLUSI abelp = sprintf('System FOM = %s', FOMPLUS) 





re(ba 

lot (3 I p) 

(ypos A _tot(shap ¢ 

dispA_tot(shape,a_cptl),"g-s', 

dispA_tot(shape,a cpt ‘b-d', 

dispA_tot(shape,a_cpt3),'r-x', 

dispA_tot(shape,a_cpt ‘m *', 

dispX_tot(shape,a_cp), ‘r--o', 

dispX_tot(shape,a_cptl),'b--s', y 

dispX_tot(shape,a_cp+2),'m-d', y 

dispX_tot(shape,a_cp+3),'c--x', ypos, .. 

dispX_tot(shape,a_cpt+4),'k--*'), grid on.. 

ae tee ives Eres %d' , ae ath 

sprintf('Re req Error = %d', error(a_cp+l)), 

sprintf('Rel Freq Error = %d', error(a_cp+2)), 

sprintf('Rel Freq Error = %d', error(a_cp+3)), 

sprintf('Rel Freq Error = %d', error(a_cp+4))); 
% legend(sprintf('Bm X, Md %d', a_cp'), sprintf('Bm X, Md %d', a_cp+l) 
% sprintf('Bm X, Md %d', a_cpt2), sprintf('Bm X, Md %d',a_cpt3),.. 
% sprintf('Bm X, Md %d', a_cp+4), sprintf('Base, Md %d', a_cp), .. 
% sprintf('Base, Md %d', a_cptl), sprintf('Base, Md %d', a_cpt2), 
% sprintf('Base, Md %d', a_cpt3), sprintf('Base, Md %d', a_cpt4) 









bar graphes of error solution using only ABC= 


ure(barp+tl0) % Tatoes 11-20 
subplot (3,2, 2*modep-1 
abc_con = int2str peen ABC(i ntervel p+modep) ) 
abc_conT = sprintf(' Cond = %s, FOM = %s', abc_con, FOMABC) 
ABC = bar(dv_cal_ABC(:,intervelptmodep),.5,'r'); hold on 
rf aankt = bar(dv_cal_ABC(:,1:5),.5,'r'); ho id’ on %trying to isolate the first 
ive modes 


base = bar(dv_cal ABC(:,1),.25,'b'); hold off%on % base first 5 modes 
%FOMabc = bar(1,0); hold off 


grid on 
Pie ea base, FOMabc], plus_conT, BASET, FOMABCl abel p), hold on 
h gr 


legend([ABC, base], abc_conT, BASET), hold on 
title(sprintf('ABC only, [ %s]', model abel p)) 


if El_Ilbls ~=0 & mass_Ibls ~=0 
% plots actual error 


stem(mass_Ibls, UL Re er age ae hold on 

stem( mass I bls, dv_mass,'k'), hold on; 

Elplot = El_Ibls+10; hold on % last half of plot 

stem Elplot, dv_El,'c','filled');hold on; 

stem El plot, dv_El,'k'); hold on 

% plots the green triangle which indicates pinned node 

%plot(barp+9,0,'g*', barp+9,0,'gh', barp+9,0,'g*',... 

% barptl1,0,'g*', barp+11,0,'gh', barpt1l1,0,'g*' ) 

plot(((BC(:,3)-1)/2),0,'g%', ((BC(:,3)-1)/2),0,'gh', ((BC(:, 3)- 
1)/2),0,'9*',.., 
((BC(:,3)-1)/2),0,'g*', ((BC(:,3)-1)/2),0,'gh', ((BC(:, 3)-1)/2),0,' g* 


elseif mass_Ibls ~=0 %&EI_Ibls =0 
% plots actual error 


stem( mass_Ibls, OS ey a ae. hold on 
stem( mass Ibls, dv_mass,'k'); hold on 
% plots the green triangle which indicates pinned node 
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%plot(barp+9,0,'g%', barpt9,0,'gh', barp+9,0,'g*') 
P ot(((BC(:,3)-1)/2),0,' 9%, ((BC(:,3)-1)/2),0,' gh’, ((BC(:, 3)- 
1)/2),0,'9*! 
else 
% plots actual error 
stem(El_Ibls, dv_El,'c','filled');hold on; 
stem El _Ibls, dv_El,'k'); hold on 
% aaa the green triangle which indicates pinned node 
F ot(barpt9,0,'g*', barp+9,0,'gh', barp+9,0,'g*') 
P ot(((BC(:,3)-1)/2),0,' 9%, ((BC(:,3)-1)/2),0,' gh’, ((BC(:, 3) 
1)/2),0,'9*! 
end % El_Ibls ... 





plus_con = int2str(cond_basePlus(intervel ptmodep)):% for legend 
} : plus_conT = sprintf('Base+ABC Cond = %s FOM = %s', plus_con, FOMPLUS); % for 
egen 
oe = eS er ae ee a ga hold on 
ase = bar(dv_cal ABC(:,1),.25,'b'); hold on % base first 5 modes 
oP OMp Ls = bar(1,0); hold off 
rid on 
Peqenatt lus, base], plus_conT, BASET), hold on 
% fegend{| plus, bas se, FOMpT us], plus_ conT, BASET, FOMPLUSI abel p), hold on 
if El_[bls ~=0 & mass_Ibls ~= 
% plots actual error 
stem( mass_Ibls, SLL a ae ely hold on 
stem( mass | bls, dv_mass,'k'); hold on 
Elplot = El_Ibls+l0; hold on % last half of plot 
stem El plot, dv_El,'c','filled'); hold on; 
stem(El plot, dv_El,'k');hold on 
% Pe the green triangle which indicates pinned node 
KE tarpnil’ Og" barpail,d, gh’, parpett: oe" gf!) 
Sots | a ee A oe MIR ne ce 
mG OF yore 
,; ((BC(:,3)-1)/2),0,'9%', ((BC(:,3)-1)/2),0,'gh', ((BC(:,3)-1)/2),0,' g* 
elseif mass_Ibls ~=0 %&EI_I bls =0 
% plots actual error 
ene eee” Gone Sehole one hold on 
stem( mass s, dv_mass,'k');ho on 
; pat the oert ae Aner ieee Paes node 
plot(barp+9,0,'g*', barp+9,0,'gh', barp+9,0,'g 
Drott({BCl:.3}-1)/2). 0) g* (CBC(:, 3)-1)/ 2) 0, ' gh! , ( (BCI 3) 
1)/2),0,'g*') 
else 
% plots actual error 
SEI Lhe: erie hela on; 
stem s, dv_El,'k'), ho on 
% plots the green triangle which indicates pinned node 
% plot barpe9o0, 9". barat. 0,. h', barp+9,0,'g*')%changed barptl 
eohanetas phot (CC BEUs 3} 11/2), 009%". (EBct e, 3-1} 72). 0, "an (ECU: 3} 
10, '9*! 


end %if El_Ibls ... 


title(sprintf('Base [1: 5] + ABC [ %s]', model abel p)) 

% title(sprintf('Base + ABC, pinned at NODE# %d', barp +1) 

a_cp= a cr +5; 
end % modep | oop 


Mogestape = mogestape + 11; % advances 

intervelp = intervelp +3; % advances to the next ABC system 
ER = ER+19; 

end % barp loop 
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hold on 


1 
1 
3 
BC Mode 1') 
2 


wtno ox won 


%ABC Mod 
% %ABC Mod 


KEK KKKKKKKKK KK KKK RKKE 


hold on 
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END plottingBARS_crs.m 


('ABC Mode 2') 
le('ABC Mode 3') 
e(' ABC Mode 5') 


1 


Oy KAKA KKK KKK KKK RK KK 


% 


V. STRAINENERGY 


HOHStrain Energy calculation of the base beam, that is neither ABC's nor 
HOHSitffnes/mass perturbations applied. 


cle 
format long 


%Given data 

% E = 106; 

% 11 = 0.016067021; 

%El ement Lenght 

L = total _length/num_el ements; 


%Making my eigenvector matrix a 42 X number of modes required 
a = zero ost 4) 2,40); 
a(3:42,:) = phi_base(:,:); 
%For loops generated to select the corresponding eigenvectors (magnitudes 
%and displacements) out of the phi_base matrix and then multiplied them 
with the corresponding hermitian Shape function second derivatives 
for i = 1:num_elements 

for k = lindof-2 
Sai = a(2*i - 1,k) 
rotleft(i,k) = a(2*i,k) 
a = a(2*i + 1,k) 
rotright(i,k) = a(2*i + 2,k) 
%"i prime" corresponds to the integral part of the strain energy equation 
%already calculated in a polynomial form, with the corresponding hermitian second 
derivative 
times its corresponding eigenvector product indicated. 
iprime(i,k) = 1/3*(6/L“J*di spleft(1,k) 42/L#rotleft(i. k)- 
6/L42*dispright(i,k)+4/L*rotright(i,k)) *3/(12/L*3*displeft(i,k) +6/L*2*rotleft(i, k)- 
12/L*3*dispright(i,k)+6/L*2*rotright(i,k))-1/3*(-6/L*2*displeft(i,k)- 
4/L*rotleft(i,k)+6/L*2*dispright(i,k)- 
2/L¥rotright(i,k))*3/(12/L*3*displeft(i,k)+6/L*2*rotleft(i,k)- 
12/L°3*dispright(i,k)+6/L*2*rotright(i,k)); 
strain(i,k) = iprime(i,k) 

end 

end 
% Assembling the strain energy equation: 

= (nominal _El/2)*strain 
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Ww. 


STRAINENERGYABC 


H%Strain Energy calculation of the base beam, with one ABC applied, 


sitffnes/mass perturbations. 


cle 
format long 


%*Gi ven data 
%E = 10e6; 
%l 1 = 0.016067021; 

L = total_length/num elements; 


a 


“This for 20 elements 


ee matrix a 42 X number of modes required 
The 42 (ndof) 


to select the corresponding eigenvectors 


base matrix and then multi p 


ie 


ral 


rresponding hermitian Shape function second deriva 


( 
li 
ti 


part of the strain ately 


orm, with the corresponding 


+6/L*2* rotleft(i,k)- 


%Making my eigen 
aa = zeros(42, 39 
Keep = [3:pinned - 1, pinned + 1:42]; 
changed if # elements changes 
aa(Keep,:) = phiaOSET(:,: 
%For loops generated 
wand displacements) out of the ph 
with the co 
for ij = 1: num_elements 

for k = lindof-3 
eae = aa(2*i - 1,k) 
rotleft(i,k) = aa(2*i,k) 
ce ee oe aa(2*i + 1,k) 
rotright(i,k) = aa(2*i + 2,k) 
%"iabcprime" corresponds to the 
already calculated in a polynom 
derivative 
times its corresponding eigenvector product indicated. 
iabcprime(i,k) = 1/3*(6/L“*2*disp 
6/L*2*dispright(i,k) +4/L*rotri ght 
12/L°3*dispright(i,k)+6/L*2*rotr 
4/L¥rotleft(i,k) +6/L*2*dispri ght ( 
2/L*rotright(i,k)) *3/(12/L*3*disp 
12/L°3*dispright(i,k)+6/L*2*rotr 


strain(i,k) = iabcprime(i, k) 
end 


end 
% Assembling the strain 
UNABC = (nominal _El/2)* ‘ 





cose 


ergy equation: 
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) 
k)); 


1 


e 


not 


quation 


is to be 


ermitian second 


~ Ir 


)- 


2*rotleft(i, 


k)- 
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